rm(list = ls())
knitr::opts_chunk$set(warnings=FALSE)
Before loading the packages install them using install()
library(tidyverse)
library(summarytools)
library(dplyr)
library(readr)
library(vegan)
library(patchwork)
library(pander)
library(ggeffects)
library(performance)
library(knitr)
library(car)
library(corrplot)
library(MASS)
library(MuMIn)
library(geosphere)
library(mgcv)
library(glmmTMB)
library(lme4)
library(factoextra)
library(fitdistrplus)
#If Issues installing packages
options("install.lock"=FALSE)
In this part the data frames were inspected in order to check the variable names and types. In some situations the data set was slightly modified for further manipulation of them.
Environmental1 <- read.csv("Data/Environmental_variables.csv")
Metadata <- read.csv("Data/Metadata_lat2018-2019.csv")
Herbivory <- read.csv("Data/Herbivory_latitudinal2018-2019.csv")
Benthic <- read.csv("Data/Percent_cover.csv")
Phlorotannins <- read.csv("Data/Phlorotannins.csv")
Predation <- read.csv("Data/Predation_latitudinal2018-2019.csv")
Size <- read.csv("Data/Size_estimations.csv")
Transects <- read.csv("Data/Transect2018-2019.csv")
Sal_data <-read.csv("Data/extracted_sal_data.csv")
temp_data <- read.csv("Data/extracted_temp_data.csv")
o2_data <- read.csv("Data/extracted_o2_data.csv")
view(dfSummary(Metadata))
#Rename temperature variables to avoid confusion with replicate measures in Environmental data
Metadata <- rename(Metadata, Temperature_site = Temperature_C)
view(dfSummary(Environmental1))
Environmental <- Environmental1 %>%
rename(D_oxygen_percent = X., D_oxygen_mgL = D_oxygen) %>% #Rename oxygen percent variable
dplyr::select(-Date,-Campaign,-Latitude) #Remove variables that are in Metadata
view(dfSummary(Benthic))
Benthic <- Benthic %>%
dplyr::select(-Date) #Remove variables that are in Metadata
view(dfSummary(Herbivory))
Herbivory <- Herbivory %>%
dplyr::select(-Date,-Campaign,-Season) %>%# Remove variables that are in Metadata
mutate(Consumed_proportion = (Initial_length_mm-Final_Length_mm)/Initial_length_mm)
It was necessary to create major categories for Bait
Herbivory <- mutate(Herbivory, Bait2 = case_when(Bait == "Lessonia_spicata" ~ "Lessonia sp.",
Bait == "Lessonia_berteroana" ~ "Lessonia sp."))
view(dfSummary(Predation))
Predation <- Predation %>%
dplyr::select(-Date,-Campaign,-Season) %>%# Remove variables that are in Metadata
mutate(Consumed_proportion = (Initial_Individuals-Remaining_individuals)/Initial_Individuals, Consumed_individuals = Initial_Individuals-Remaining_individuals) #Add Consumed_proportion column
Create major categories for Bait (Squidpops vs Crabs) and Assay time (2 hours and 24 hours)
Predation <- mutate(Predation, Bait2 = case_when(Bait == "Squidpops" ~ "Squidpops",
Bait != "Squidpops" ~ "Crabs"), Assay_time2 = case_when(Assay_time == "24" ~ "24",
Bait != "24" ~ "2"))
view(dfSummary(Phlorotannins))
Phlorotannins <- Phlorotannins %>%
dplyr::select(-Date,-Latitude) # Remove variables that are in metadata
view(dfSummary(Size))
Size <- Size %>%
dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata
view(dfSummary(Transects))
Transects <- Transects %>%
dplyr::select(-Date,-Campaign) # Remove variables that are in Metadata
In this section data sets are manipulated for analyses. New data sets with means and combined with metadata and consumption are created.
Environmental_metadata <- Environmental %>%
left_join(Metadata, by = "Site")
Environmental_2 <- Environmental_metadata %>%
group_by(Site, Latitude) %>%
dplyr::select(-D_oxygen_percent, -Conductivity_Ms, -C_f, -C_i, -H_Visibility) #Keep only variables that will be used in the analysis
Environmental_2$Day_length <- daylength(Environmental_2$Latitude, Environmental_2$Date)
Environmental_sum_long <- Environmental_2 %>%
dplyr::select(Site, Latitude, Temperature, D_oxygen_mgL, Salinity_ppt, V_Visibility, Day_length) %>%
gather(Env_Variable, Value, Temperature:Day_length) %>%
group_by(Site, Latitude, Env_Variable)%>%
summarise(Mean_value = mean(Value, na.rm=T),n = n(), sd = sd(Value,na.rm=T), se = sd/sqrt(n))
Environmental_sum <- Environmental_sum_long %>%
dplyr::select(Site, Env_Variable, Mean_value) %>%
pivot_wider(names_from = Env_Variable, values_from = Mean_value,
values_fill = 0)
Then, add environmental data to consumption
Predation vs environmental data sets
Predation_env <- Predation %>%
left_join(Environmental_sum, by = "Site")
Herbivory_env <- Herbivory %>%
left_join(Environmental_sum, by = "Site")
Benthic_metadata <- Benthic %>%
left_join(Metadata, by = "Site")
Benthic_metadata_wide <- Benthic_metadata %>%
group_by(Site, Latitude, Type) %>%
pivot_wider(names_from = Type, values_from = Coverage,
values_fill = 0)
Benthic_mean <- Benthic_metadata %>%
group_by(Site, Latitude, Type)%>%
summarise(Mean_cover = mean(Coverage),n = n(), sd = sd(Coverage), se = sd/sqrt(n))
Benthic_sum_wide <- Benthic %>%
group_by(Site, Type) %>%
summarise(cover_mean = mean(Coverage)) %>%
pivot_wider(names_from = Type, values_from = cover_mean,
values_fill = 0)
Benthic_sum_wide_var <- dplyr::select(Benthic_sum_wide, Site, Kelp, Bare_rock, Understorey_algae)
Then, add benthic data to consumption
Predation_env + Benthic summarised data
Predation_2 <- Predation_env %>%
left_join(Benthic_sum_wide_var, by = "Site")
Herbivory_2 <- Herbivory_env %>%
left_join(Benthic_sum_wide_var, by = "Site")
Size_metadata <- Size %>%
left_join(Metadata, by = c("Site"))
Size_mean_sp<- Size_metadata %>%
group_by(Site, Latitude, Species, Type)%>%
summarise(mean_Size_sp = mean(Size_cm))
Size_mean <- Size_mean_sp %>%
group_by(Site, Latitude, Type)%>%
summarise(mean_Size = mean(mean_Size_sp),n = n(), sd = sd(mean_Size_sp), se = sd/sqrt(n))
Size_mean_wide <- Size_mean_sp %>%
filter(Type == "Fish") %>%
group_by(Site, Type)%>%
summarise(mean_Size = mean(mean_Size_sp)) %>%
pivot_wider(names_from = Type, values_from = mean_Size,
values_fill = 0) %>%
rename(Fish_size = Fish)
Add fish length data to consumption
Predation_env + fish size data
Predation_3 <- Predation_2 %>%
left_join(Size_mean_wide, by = "Site")
#Join transects with metadata
Transect_metadata <- Transects %>%
left_join(Metadata, by = c("Site"))
#Summary of species from transects
Transect_sum <- Transect_metadata %>%
filter(Method == "Transect") %>% #Select only data taken with visual transect method
group_by(Site, Latitude, Census, Species) %>%
summarise(Total_abundance = sum(Individuals))
First, change data frame to the wide format
Transect_wide <- Transect_sum %>%
ungroup() %>%
pivot_wider(names_from = Species, values_from = Total_abundance,
values_fill = 0)
Convert back to long format with Species and Abundance as columns. This will add rows with zeroes for all the species missing in each video
Transect_long <- Transect_wide %>%
gather(Species, Abundance, Anisotremus_scapularis:Labrisomus_philippii)
#Obtain mean abundance per site and method
Transect_mean <- Transect_long %>%
group_by(Site, Latitude, Species) %>%
summarise(Mean_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))
Follow the same steps as abundances by type
TransectGroup_sum <- Transect_metadata %>%
group_by(Site, Latitude, Census, Group) %>%
summarise(Total_abundance = sum(Individuals))
TransectGroup_wide <- TransectGroup_sum %>%
ungroup() %>%
pivot_wider(names_from = Group, values_from = Total_abundance,
values_fill = 0)
TransectGroup_wide <- TransectGroup_wide%>%
mutate(Herbivores = Herbivore_fish + Herbivore_gastropod + Urchin + Omnivorous_crab, Invertebrate_predators = Predatory_gastropod+ Snail + Predatory_crab + Sea_star + Omnivorous_crab + Schrimp)
TransectGroup_long <- TransectGroup_wide %>%
gather(Group, Abundance, Herbivore_fish:Invertebrate_predators)
TransectGroup_mean <- TransectGroup_long %>%
group_by(Site, Latitude, Group) %>%
summarise(Total_abundance = mean(Abundance),n = n(), sd = sd(Abundance), se = sd/sqrt(n))
First, add a column of presence when abundance is higher than 0.
Transect_long <- Transect_long %>%
mutate(Presence = if_else(Abundance > 0,1,0))
Calculate sum of all species in the same site and census.
Richness_transect <- Transect_long %>%
group_by(Site, Latitude, Census) %>%
summarise(Richness = sum(Presence))
#Calculate mean by site
Richness_mean <- Richness_transect %>%
group_by(Site, Latitude)%>%
summarise(Mean_SR = mean(Richness),n = n(), sd = sd(Richness), se = sd/sqrt(n))
Calculate means and change it to the wide format
#Per Group
Group_abundances_mean <- TransectGroup_long%>%
filter(Group == "Invertebrate_predators" | Group == "Predatory_fish"| Group == "Herbivores")%>%
group_by(Site, Group)%>%
summarise(Mean_Abundance = mean(Abundance)) %>%
pivot_wider(names_from = Group, values_from = Mean_Abundance,
values_fill = 0)
Then, add abundances data to consumption
#Select only the needed columns
Group_abundances_mean_pred <- dplyr::select(Group_abundances_mean, Site, Invertebrate_predators, Predatory_fish)
Predation_4 <- Predation_3 %>%
left_join(Group_abundances_mean_pred, by = "Site")
#Select only the needed columns
Group_abundances_mean_herb <- dplyr::select(Group_abundances_mean, Site, Herbivores)
Herbivory_3 <- Herbivory_2 %>%
left_join(Group_abundances_mean_herb, by = "Site")
Calculate means and change it to the wide format
Richness_meanpred <- Richness_transect %>%
group_by(Site)%>%
summarise(Mean_SR = mean(Richness))
Then, add abundances data to consumption
Predation_5 <- Predation_4 %>%
left_join(Richness_meanpred, by = "Site")
Herbivory_4 <- Herbivory_3 %>%
left_join(Richness_meanpred, by = "Site")
Summarise abundances of species per site and standardise
Transect_sum <- Transect_long %>%
group_by(Site, Latitude, Species) %>%
summarise(Total_abundance = mean(Abundance))
## `summarise()` has grouped output by 'Site', 'Latitude'. You can override using
## the `.groups` argument.
Transect_wide2 <- Transect_sum %>%
ungroup() %>%
pivot_wider(names_from = Species, values_from = Total_abundance, values_fill = 0)
Transect_val <- Transect_wide2 %>%
dplyr::select(- Site, -Latitude)
Transect_stand <- decostand(Transect_val, method = "total")
Calculate Bray-Curtis dissimilariry
Transect_comp <- vegdist(Transect_stand, "bray")
Run Principal Coordinate Analysis (PcoA) to obtain the component explaining the highers variation
pcoa <- cmdscale(Transect_comp, eig = TRUE, add = TRUE)
pcoa$eig
## [1] 1.740530e+00 1.329554e+00 8.438524e-01 6.203963e-01 5.740214e-01
## [6] 3.517972e-01 1.671264e-01 7.777519e-02 5.878784e-02 3.878297e-02
## [11] -7.537429e-17 -2.615278e-16
(pcoa$eig/sum(pcoa$eig))*100 #proportion of variation
## [1] 2.999557e+01 2.291298e+01 1.454260e+01 1.069165e+01 9.892446e+00
## [6] 6.062727e+00 2.880187e+00 1.340345e+00 1.013125e+00 6.683695e-01
## [11] -1.298969e-15 -4.507061e-15
#convert pcoa results into data frame that can be plotted
pcoa_df <- data.frame(pcoa$points)
colnames(pcoa_df) <- c("PCo1", "PCo2")
pcoa_df <- pcoa_df %>%
dplyr::select(-c(PCo2))
pcoa_df$Site <- factor(Transect_wide2$Site)
pcoa_df$Latitude <- Transect_wide2$Latitude
PCo 1 represents 30% of the variation community composition
pcoa_df1<- dplyr::select(pcoa_df, -Latitude)
Predation_6 <- Predation_5 %>%
left_join(pcoa_df1, by = "Site")
Herbivory_5 <- Herbivory_4 %>%
left_join(pcoa_df1, by = "Site")
Phlorotannins_metadata <- Phlorotannins %>%
left_join(Metadata, by = c("Site"))
Phlorotannins_mean <- Phlorotannins_metadata %>%
group_by(Site, Latitude) %>%
summarise(phloro_mean = mean(Phlorotannins), n = n(), sd = sd(Phlorotannins), se = sd/sqrt(n))
Phlorotannins_sum <- Phlorotannins %>%
group_by(Site) %>%
summarise(phloro_mean = mean(Phlorotannins))
Herbivory_6 <- Herbivory_5 %>%
left_join(Phlorotannins_sum, by = "Site")
Latitudinal patterns were analysed for the recorded explanatory variables using linear models.
Set general theme for plots
Theme1 <- theme_bw() +
theme(panel.grid.major = element_line(colour = "white"),
panel.grid.minor = element_line(colour = "white"), axis.text.y = element_text(size = 9), axis.title.y = element_text(size = 9), axis.text.x = element_text(size = 9), axis.title.x = element_text(size = 9))
For these analyses we used the approach and code of model selection for environmental and spatial gradients provided by Anderson et al. (2022)
a. Analyses
For future plotting
Environmental_2$Latitude <- Environmental_2$Latitude*-1
#change sign of latitude
# First, examine a simple scatter plot (always a good idea!)
plot(Environmental_2$Latitude, Environmental_2$Temperature, xlab = "Latitude", ylab = "Temperature", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Environmental_2$Latitude, Environmental_2$Temperature, pch = 21, col = mygrey, bg = mygrey)
descdist(Environmental_2$Temperature, discrete = FALSE)
## summary statistics
## ------
## min: 10.2 max: 18.7
## median: 16
## mean: 15.32917
## estimated sd: 2.365349
## estimated skewness: -0.7647711
## estimated kurtosis: 2.651633
normal_dist <- fitdist(Environmental_2$Temperature, "weibull")
plot(normal_dist)
T_mod1 <- glmer(Temperature ~ Latitude+ (1|Site), data= Environmental_2, family = Gamma(link = "log"))
T_mod2 <- glmer(Temperature ~poly(Latitude,2)+ (1|Site), data= Environmental_2, family = Gamma(link = "log"))
T_mod3 <- glmer(Temperature ~poly(Latitude,3)+ (1|Site), data= Environmental_2, family = Gamma(link = "log"))
AICc(T_mod1, T_mod3, T_mod2)
## df AICc
## T_mod1 4 80.68285
## T_mod3 6 84.78450
## T_mod2 5 82.55903
r2(T_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.930
## Marginal R2: 0.713
r2(T_mod2)
## # R2 for Mixed Models
##
## Conditional R2: 0.935
## Marginal R2: 0.744
r2(T_mod3)
## # R2 for Mixed Models
##
## Conditional R2: 0.937
## Marginal R2: 0.761
T_mod1 is the best fitted model
res <- resid(T_mod1)
#produce residual vs. fitted plot
plot(fitted(T_mod1), res)
#add a horizontal line at 0
abline(0,0)
#create Q-Q plot for residuals
qqnorm(res)
#add a straight diagonal line to the plot
qqline(res)
summary(T_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: Temperature ~ Latitude + (1 | Site)
## Data: Environmental_2
##
## AIC BIC logLik deviance df.resid
## 79.8 87.2 -35.9 71.8 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.34203 -0.41973 0.05954 0.42096 1.54302
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.004483 0.06696
## Residual 0.001452 0.03811
## Number of obs: 48, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 3.23945 0.20465 15.829 < 2e-16 ***
## Latitude -0.01755 0.00669 -2.623 0.00871 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Latitude -0.975
r2(T_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.930
## Marginal R2: 0.713
b. Plot predicted relationship
predictions <- ggpredict(T_mod1, terms ="Latitude")
# Plot the GAM model predictions
Temperature_plot <- ggplot() +
geom_jitter(data = Environmental_2, aes(x= Latitude, y=Temperature), width = 0.08, height = 0.03, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.high, ymax = conf.low), alpha = 0.2) +
geom_line(data = predictions, aes(x = x, y =predicted), linewidth = 0.6) +
Theme1 +
labs( y="Temperature (°C)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank())
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(Environmental_2$Latitude, Environmental_2$Salinity_ppt, xlab = "Latitude", ylab = "Salinity", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Environmental_2$Latitude, Environmental_2$Salinity_ppt, pch = 21, col = mygrey, bg = mygrey)
descdist(Environmental_2$Salinity_ppt, discrete = FALSE)
## summary statistics
## ------
## min: 32.8 max: 34.4
## median: 33.8
## mean: 33.73125
## estimated sd: 0.3882071
## estimated skewness: -0.9322086
## estimated kurtosis: 3.328625
normal_dist <- fitdist(Environmental_2$Salinity_ppt, "weibull")
plot(normal_dist)
S_mod1 <- glmer(Salinity_ppt ~ Latitude+ (1|Site), data= Environmental_2, family = Gamma(link = "log"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0215438 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
S_mod2 <- glmer(Salinity_ppt ~poly(Latitude,2)+ (1|Site), data= Environmental_2, family = Gamma(link = "log"))
S_mod3 <- glmer(Salinity_ppt ~poly(Latitude,3)+ (1|Site), data= Environmental_2, family = Gamma(link = "log"))
check_model(S_mod1)
AICc(S_mod1, S_mod3, S_mod2)
## df AICc
## S_mod1 4 20.65664
## S_mod3 6 24.12187
## S_mod2 5 21.50425
BIC(S_mod1, S_mod3, S_mod2)
## df BIC
## S_mod1 4 27.21121
## S_mod3 6 33.30029
## S_mod2 5 29.43168
summary(S_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: Salinity_ppt ~ Latitude + (1 | Site)
## Data: Environmental_2
##
## AIC BIC logLik deviance df.resid
## 19.7 27.2 -5.9 11.7 44
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8003 -0.5564 0.1690 0.5914 1.4235
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 2.189e-05 0.004679
## Residual 5.926e-05 0.007698
## Number of obs: 48, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 3.5487577 0.0083718 423.896 < 2e-16 ***
## Latitude -0.0010191 0.0002757 -3.696 0.000219 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Latitude -0.953
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.0215438 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?
r2(S_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.547
## Marginal R2: 0.379
b. Plot predicted relationship
predictions <- ggpredict(S_mod1, terms ="Latitude [all]")
# Plot the GAM model predictions
Salinity_plot <- ggplot() +
geom_jitter(data = Environmental_2, aes(x= Latitude, y=Salinity_ppt), width = 0.08, height = 0.05, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
geom_line(data =predictions, aes(x = x, y = predicted), linewidth = 0.6) +
Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank())
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, xlab = "Latitude", ylab = "DO", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Environmental_2$Latitude, Environmental_2$D_oxygen_mgL, pch = 21, col = mygrey, bg = mygrey)
descdist(Environmental_2$D_oxygen_mgL, discrete = FALSE)
## summary statistics
## ------
## min: 1.62 max: 12.15
## median: 6.445
## mean: 6.381042
## estimated sd: 2.626742
## estimated skewness: 0.3329725
## estimated kurtosis: 2.893175
normal_dist <- fitdist(Environmental_2$D_oxygen_mgL, "norm")
plot(normal_dist)
Ox_mod1 <- glmer(D_oxygen_mgL ~ Latitude+ (1|Site), data= Environmental_2, family = gaussian)
## Warning in glmer(D_oxygen_mgL ~ Latitude + (1 | Site), data = Environmental_2,
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
Ox_mod2 <- glmer(D_oxygen_mgL ~poly(Latitude,2)+ (1|Site), data= Environmental_2, family = gaussian)
## Warning in glmer(D_oxygen_mgL ~ poly(Latitude, 2) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
Ox_mod3 <- glmer(D_oxygen_mgL ~poly(Latitude,3)+ (1|Site), data= Environmental_2, family = gaussian)
## Warning in glmer(D_oxygen_mgL ~ poly(Latitude, 3) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
AICc(Ox_mod1, Ox_mod2, Ox_mod3)
## df AICc
## Ox_mod1 4 221.6313
## Ox_mod2 5 209.4476
## Ox_mod3 6 207.2213
BIC(Ox_mod1, Ox_mod2, Ox_mod3)
## df BIC
## Ox_mod1 4 228.1858
## Ox_mod2 5 217.3750
## Ox_mod3 6 216.3998
Improvement cubic model shows the best fit
Two options may be possible based on AICc. Very low df in Ox_mod 3
summary(Ox_mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: D_oxygen_mgL ~ poly(Latitude, 3) + (1 | Site)
## Data: Environmental_2
##
## REML criterion at convergence: 193.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37338 -0.33408 0.05551 0.26713 3.04460
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 4.611 2.147
## Residual 3.029 1.741
## Number of obs: 48, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 6.3810 0.6688 9.540
## poly(Latitude, 3)1 0.6637 4.6339 0.143
## poly(Latitude, 3)2 -6.5554 4.6339 -1.415
## poly(Latitude, 3)3 -0.1882 4.6339 -0.041
##
## Correlation of Fixed Effects:
## (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000
## ply(Ltt,3)2 0.000 0.000
## ply(Ltt,3)3 0.000 0.000 0.000
#Plot the GAM model predictions
Oxygen_plot <- ggplot() +
geom_jitter(data = Environmental_2, aes(x= Latitude, y=D_oxygen_mgL), width = 0.08, height = 0.07, size = 1, alpha = 0.6) +
Theme1 +
labs( y="Dissolved oxygen (mg"~L^-1*")", x="Latitude (°S)") +
theme(axis.title.x=element_text(size = 10),
axis.text.x=element_text(size = 10), axis.text.y = element_text(size = 10), axis.title.y = element_text(size = 8))
par(mfrow = c(2, 2))
check_model(Ox_mod3)
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(Environmental_2$Latitude, Environmental_2$V_Visibility, xlab = "Latitude", ylab = "Visibility", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Environmental_2$Latitude, Environmental_2$V_Visibility, pch = 21, col = mygrey, bg = mygrey)
descdist(Environmental_2$V_Visibility, discrete = FALSE)
## summary statistics
## ------
## min: 3 max: 10
## median: 5.75
## mean: 5.854167
## estimated sd: 1.994563
## estimated skewness: 0.4867126
## estimated kurtosis: 2.512658
normal_dist <- fitdist(Environmental_2$V_Visibility, "norm")
plot(normal_dist)
V_mod1 <- glmer(V_Visibility ~ Latitude+ (1|Site), data= Environmental_2, family = gaussian())
## Warning in glmer(V_Visibility ~ Latitude + (1 | Site), data = Environmental_2,
## : calling glmer() with family=gaussian (identity link) as a shortcut to lmer()
## is deprecated; please call lmer() directly
V_mod2 <- glmer(V_Visibility ~poly(Latitude,2)+ (1|Site), data= Environmental_2, family = gaussian())
## Warning in glmer(V_Visibility ~ poly(Latitude, 2) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
V_mod3 <- glmer(V_Visibility ~poly(Latitude,3)+ (1|Site), data= Environmental_2, family = gaussian())
## Warning in glmer(V_Visibility ~ poly(Latitude, 3) + (1 | Site), data =
## Environmental_2, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
check_model(V_mod1)
AICc(V_mod1, V_mod2, V_mod3)
## df AICc
## V_mod1 4 171.7591
## V_mod2 5 162.0855
## V_mod3 6 158.9580
BIC(V_mod1, V_mod2, V_mod3)
## df BIC
## V_mod1 4 178.3137
## V_mod2 5 170.0129
## V_mod3 6 168.1365
V_mod3 is the best fitted
summary(V_mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: V_Visibility ~ poly(Latitude, 3) + (1 | Site)
## Data: Environmental_2
##
## REML criterion at convergence: 144.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.83174 -0.54808 0.07988 0.44641 2.45295
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 3.2862 1.8128
## Residual 0.8715 0.9336
## Number of obs: 48, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.8542 0.5404 10.834
## poly(Latitude, 3)1 5.0601 3.7438 1.352
## poly(Latitude, 3)2 -0.1325 3.7438 -0.035
## poly(Latitude, 3)3 -4.2251 3.7438 -1.129
##
## Correlation of Fixed Effects:
## (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000
## ply(Ltt,3)2 0.000 0.000
## ply(Ltt,3)3 0.000 0.000 0.000
# Plot the GAM model predictions
Visibility_plot <- ggplot() +
geom_jitter(data = Environmental_2, aes(x= Latitude, y=V_Visibility), width = 0.08, height = 0.08, size = 1, alpha = 0.6) +
Theme1 +
labs( y="Visibility (m)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank())
a. Analyses
pcoa_df$Latitude <- pcoa_df$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
plot(pcoa_df$Latitude, pcoa_df$PCo1, xlab = "Latitude", ylab = "PCoA1", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(pcoa_df$Latitude, pcoa_df$PCo1, pch = 21, col = mygrey, bg = mygrey)
descdist(pcoa_df$PCo1, discrete = FALSE)
## summary statistics
## ------
## min: -0.3828236 max: 0.6887275
## median: -0.1217363
## mean: 5.247539e-17
## estimated sd: 0.3977814
## estimated skewness: 0.9122292
## estimated kurtosis: 2.331104
normal_dist <- fitdist(log(pcoa_df$PCo1+1), "norm")
plot(normal_dist)
PCO_mod1 <- glm(log(PCo1 +1) ~ Latitude, data= pcoa_df, family = gaussian())
PCO_mod2 <- glm(log(PCo1 +1) ~poly(Latitude,2), data= pcoa_df, family = gaussian())
PCO_mod3 <- glm(log(PCo1 +1) ~poly(Latitude,3), data= pcoa_df, family = gaussian())
check_model(PCO_mod2)
check_model(PCO_mod1)
AICc(PCO_mod1,PCO_mod2, PCO_mod3)
## df AICc
## PCO_mod1 3 14.72754
## PCO_mod2 4 17.28033
## PCO_mod3 5 23.21428
summary(PCO_mod1)
##
## Call:
## glm(formula = log(PCo1 + 1) ~ Latitude, family = gaussian(),
## data = pcoa_df)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.85098 0.43453 -1.958 0.0786 .
## Latitude 0.02632 0.01421 1.853 0.0936 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.1132384)
##
## Null deviance: 1.5211 on 11 degrees of freedom
## Residual deviance: 1.1324 on 10 degrees of freedom
## AIC: 11.728
##
## Number of Fisher Scoring iterations: 2
par(mfrow = c(2, 2))
gam.check(PCO_mod1)
##
## 'gamm' based fit - care required with interpretation.
## Checks based on working residuals may be misleading.
No significant effect
# Plot the GAM model predictions
PCOA_plot <- ggplot() +
geom_jitter(data = pcoa_df, aes(x= Latitude, y=PCo1), width = 0, height = 0, size = 1, alpha = 0.6) +
Theme1 +
labs( y="log(PCoA1 + 1)", x="Latitude (°S)")
a. Analyses
For further plotting
Benthic_metadata_wide$Latitude <- Benthic_metadata_wide$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, xlab = "Latitude", ylab = "Kelp", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Kelp, pch = 21, col = mygrey, bg = mygrey)
Kelp_mod1 <- glmmTMB(Kelp ~ Latitude+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
Kelp_mod2 <- glmmTMB(Kelp ~poly(Latitude,2)+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
Kelp_mod3 <- glmmTMB(Kelp ~poly(Latitude,3)+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
check_model(Kelp_mod3)
## `check_outliers()` does not yet support models of class `glmmTMB`.
AICc(Kelp_mod2, Kelp_mod1, Kelp_mod3)
## df AICc
## Kelp_mod2 6 470.8080
## Kelp_mod1 5 468.5267
## Kelp_mod3 7 473.2009
BIC(Kelp_mod2, Kelp_mod1, Kelp_mod3)
## df BIC
## Kelp_mod2 6 483.7652
## Kelp_mod1 5 479.4769
## Kelp_mod3 7 488.0979
summary(Kelp_mod1)
## Family: tweedie ( log )
## Formula: Kelp ~ Latitude + (1 | Site)
## Data: Benthic_metadata_wide
##
## AIC BIC logLik deviance df.resid
## 467.7 479.5 -228.8 457.7 73
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 1.159 1.076
## Number of obs: 78, groups: Site, 12
##
## Dispersion parameter for tweedie family (): 6.27
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.42602 1.70405 3.771 0.000163 ***
## Latitude -0.15074 0.06016 -2.506 0.012217 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(Kelp_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.800
## Marginal R2: 0.371
b. Plot predicted relationship
predictions <- ggpredict(Kelp_mod1, terms ="Latitude [all]")
# Plot the GAM model predictions
Kelp_plot <- ggplot() +
geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Kelp), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) +
Theme1 +
labs( y="Kelp cover (%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y=element_text(8))
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, xlab = "Latitude", ylab = "Bare_rock", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Bare_rock, pch = 21, col = mygrey, bg = mygrey)
BR_mod1 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma= 1.4, data = Benthic_metadata_wide, family = gaussian())
BR_mod2 <- mgcv::gam(Bare_rock ~ s(Latitude), gamma=1.4, data = Benthic_metadata_wide, family = tw())
descdist(Benthic_metadata_wide$Bare_rock, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 100
## median: 23.85013
## mean: 36.00036
## estimated sd: 32.50373
## estimated skewness: 0.5844581
## estimated kurtosis: 1.884422
normal_dist <- fitdist(Benthic_metadata_wide$Bare_rock, "exp")
plot(normal_dist)
BR_mod1 <- glmmTMB(Bare_rock ~ Latitude+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
BR_mod2 <- glmmTMB(Bare_rock ~poly(Latitude,2)+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
BR_mod3 <- glmmTMB(Bare_rock ~poly(Latitude,3)+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
check_model(BR_mod2)
## `check_outliers()` does not yet support models of class `glmmTMB`.
AICc(BR_mod1, BR_mod2, BR_mod3)
## df AICc
## BR_mod1 5 634.5401
## BR_mod2 6 634.5044
## BR_mod3 7 636.9098
BIC(BR_mod1, BR_mod2, BR_mod3)
## df BIC
## BR_mod1 5 645.4903
## BR_mod2 6 647.4616
## BR_mod3 7 651.8068
summary(BR_mod2)
## Family: tweedie ( log )
## Formula: Bare_rock ~ poly(Latitude, 2) + (1 | Site)
## Data: Benthic_metadata_wide
##
## AIC BIC logLik deviance df.resid
## 633.3 647.5 -310.7 621.3 72
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.8647 0.9299
## Number of obs: 78, groups: Site, 12
##
## Dispersion parameter for tweedie family (): 4.1
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.126 0.283 11.046 <2e-16 ***
## poly(Latitude, 2)1 2.890 2.449 1.180 0.238
## poly(Latitude, 2)2 -3.909 2.433 -1.606 0.108
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(BR_mod2)
## # R2 for Mixed Models
##
## Conditional R2: 0.859
## Marginal R2: 0.225
R2 of mod2 is higher
# Plot the GAM model predictions
Bare_plot <- ggplot() +
geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Bare_rock), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
Theme1 +
labs( y="Bare rock cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y=element_text(8))
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Understorey_algae, xlab = "Latitude", ylab = "Understorey_algae", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Benthic_metadata_wide$Latitude, Benthic_metadata_wide$Understorey_algae, pch = 21, col = mygrey, bg = mygrey)
descdist(Benthic_metadata_wide$Understorey_algae, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 82.92683
## median: 0
## mean: 10.92887
## estimated sd: 19.35412
## estimated skewness: 2.110129
## estimated kurtosis: 7.04205
normal_dist <- fitdist(Benthic_metadata_wide$Understorey_algae, "exp")
plot(normal_dist)
UA_mod1 <- glmmTMB(Understorey_algae ~ Latitude+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
UA_mod2 <- glmmTMB(Understorey_algae ~poly(Latitude,2)+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
UA_mod3 <- glmmTMB(Understorey_algae ~poly(Latitude,3)+ (1|Site), data= Benthic_metadata_wide, family = tweedie())
check_model(UA_mod2)
## `check_outliers()` does not yet support models of class `glmmTMB`.
AICc(UA_mod1, UA_mod2, UA_mod3)
## df AICc
## UA_mod1 5 342.5629
## UA_mod2 6 344.5077
## UA_mod3 7 342.1709
BIC(UA_mod1, UA_mod2, UA_mod3)
## df BIC
## UA_mod1 5 353.5131
## UA_mod2 6 357.4649
## UA_mod3 7 357.0679
summary(UA_mod1)
## Family: tweedie ( log )
## Formula: Understorey_algae ~ Latitude + (1 | Site)
## Data: Benthic_metadata_wide
##
## AIC BIC logLik deviance df.resid
## 341.7 353.5 -165.9 331.7 73
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.819 2.412
## Number of obs: 78, groups: Site, 12
##
## Dispersion parameter for tweedie family (): 6.61
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.75831 3.33416 0.527 0.598
## Latitude -0.04022 0.11085 -0.363 0.717
r2(UA_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.934
## Marginal R2: 0.011
# Plot the GAM model predictions
UA_plot <- ggplot() +
geom_jitter(data = Benthic_metadata_wide, aes(x= Latitude, y=Understorey_algae), width = 0.08, height = 0.11, size = 1, alpha = 0.6) +
Theme1 +
labs( y="UA cover(%)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y=element_text(8))
TransectGroup_wide$Latitude <- TransectGroup_wide$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
plot(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, xlab = "Latitude", ylab = "PF", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(TransectGroup_wide$Latitude, TransectGroup_wide$Predatory_fish, pch = 21, col = mygrey, bg = mygrey)
descdist(TransectGroup_wide$Predatory_fish, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 333
## median: 29
## mean: 51.61538
## estimated sd: 75.31432
## estimated skewness: 2.720417
## estimated kurtosis: 10.84418
normal_dist <- fitdist(TransectGroup_wide$Predatory_fish, "exp")
plot(normal_dist)
PF_mod1 <- glmer.nb(Predatory_fish ~ Latitude+ (1|Site), data= TransectGroup_wide)
PF_mod2 <- glmer.nb(Predatory_fish ~poly(Latitude,2)+ (1|Site), data= TransectGroup_wide)
PF_mod3 <- glmer.nb(Predatory_fish ~poly(Latitude,3)+ (1|Site), data= TransectGroup_wide)
## boundary (singular) fit: see help('isSingular')
check_model(PF_mod3)
AICc(PF_mod2, PF_mod1, PF_mod3)
## df AICc
## PF_mod2 5 257.1651
## PF_mod1 4 257.3351
## PF_mod3 6 249.9318
BIC(PF_mod2, PF_mod1, PF_mod3)
## df BIC
## PF_mod2 5 260.4556
## PF_mod1 4 260.4627
## PF_mod3 6 253.0594
mod 3 is singular, so mod 1 was chosen
summary(PF_mod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1.2985) ( log )
## Formula: Predatory_fish ~ poly(Latitude, 2) + (1 | Site)
## Data: TransectGroup_wide
##
## AIC BIC logLik deviance df.resid
## 254.2 260.5 -122.1 244.2 21
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0949 -0.6533 -0.2637 0.1770 2.9381
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4711 0.6864
## Number of obs: 26, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3709 0.3024 11.147 < 2e-16 ***
## poly(Latitude, 2)1 -4.8871 1.5132 -3.230 0.00124 **
## poly(Latitude, 2)2 -2.7915 1.5243 -1.831 0.06705 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(L,2)1
## ply(Ltt,2)1 0.071
## ply(Ltt,2)2 0.159 0.175
r2(PF_mod2)
## # R2 for Mixed Models
##
## Conditional R2: 0.749
## Marginal R2: 0.546
b. Plot predicted relationship
predictions <- ggpredict(PF_mod2, terms ="Latitude [all]")
# Plot the GAM model predictions
PredFish_plot <- ggplot() +
geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Predatory_fish), width = 0.12, height = 0.5, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) +
Theme1 +
labs( y="Predatory Fish", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y=element_text(8))
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, xlab = "Latitude", ylab = "Herbivores", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(TransectGroup_wide$Latitude, TransectGroup_wide$Herbivores, pch = 21, col = mygrey, bg = mygrey)
descdist(TransectGroup_wide$Herbivores, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 1912
## median: 155.5
## mean: 389.4231
## estimated sd: 547.3004
## estimated skewness: 1.702098
## estimated kurtosis: 5.106245
normal_dist <- fitdist(TransectGroup_wide$Herbivores, "exp")
plot(normal_dist)
Herb_mod1 <- glmer.nb(Herbivores ~ Latitude+ (1|Site), data= TransectGroup_wide)
Herb_mod2 <- glmer.nb(Herbivores ~poly(Latitude,2)+ (1|Site), data= TransectGroup_wide)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0495445 (tol = 0.002, component 1)
Herb_mod3 <- glmer.nb(Herbivores ~poly(Latitude,3)+ (1|Site), data= TransectGroup_wide)
check_model(Herb_mod1)
AICc(Herb_mod1, Herb_mod2, Herb_mod3)
## df AICc
## Herb_mod1 4 326.4081
## Herb_mod2 5 326.9362
## Herb_mod3 6 327.9244
BIC(Herb_mod1, Herb_mod2, Herb_mod3)
## df BIC
## Herb_mod1 4 329.5358
## Herb_mod2 5 330.2266
## Herb_mod3 6 331.0519
summary(Herb_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1.5227) ( log )
## Formula: Herbivores ~ Latitude + (1 | Site)
## Data: TransectGroup_wide
##
## AIC BIC logLik deviance df.resid
## 324.5 329.5 -158.3 316.5 22
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.2022 -0.4869 -0.0675 0.3265 1.3791
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.748 2.398
## Number of obs: 26, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 10.1339 3.3117 3.060 0.00221 **
## Latitude -0.2029 0.1101 -1.844 0.06523 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Latitude -0.976
r2(Herb_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.939
## Marginal R2: 0.250
# Plot the GAM model predictions
Herb_plot <- ggplot() +
geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Herbivores), width = 0.2, height = 0.5, size = 1, alpha = 0.6) +
Theme1 +
labs( y="Herbivores", x="Latitude (°S)")
# First, examine a simple scatter plot (always a good idea!)
plot(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, xlab = "Latitude", ylab = "Invertebrate pred", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(TransectGroup_wide$Latitude, TransectGroup_wide$Invertebrate_predators, pch = 21, col = mygrey, bg = mygrey)
descdist(TransectGroup_wide$Invertebrate_predators, discrete = FALSE)
## summary statistics
## ------
## min: 8 max: 753
## median: 174.5
## mean: 242.3077
## estimated sd: 208.9332
## estimated skewness: 1.010871
## estimated kurtosis: 3.417233
normal_dist <- fitdist(TransectGroup_wide$Invertebrate_predators, "exp")
plot(normal_dist)
Inv_mod1 <- glmer.nb(Invertebrate_predators ~ Latitude+ (1|Site), data= TransectGroup_wide)
Inv_mod2 <- glmer.nb(Invertebrate_predators ~poly(Latitude,2)+ (1|Site), data= TransectGroup_wide)
Inv_mod3 <- glmer.nb(Invertebrate_predators ~poly(Latitude,3)+ (1|Site), data= TransectGroup_wide)
check_model(Inv_mod1)
AICc(Inv_mod2, Inv_mod1, Inv_mod3)
## df AICc
## Inv_mod2 5 337.9285
## Inv_mod1 4 337.7932
## Inv_mod3 6 337.1040
BIC(Inv_mod2, Inv_mod1, Inv_mod3)
## df BIC
## Inv_mod2 5 341.2190
## Inv_mod1 4 340.9208
## Inv_mod3 6 340.2315
summary(Inv_mod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(3.6522) ( log )
## Formula: Invertebrate_predators ~ poly(Latitude, 3) + (1 | Site)
## Data: TransectGroup_wide
##
## AIC BIC logLik deviance df.resid
## 332.7 340.2 -160.3 320.7 20
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.48939 -0.60016 -0.01322 0.64400 1.52814
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4727 0.6875
## Number of obs: 26, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.0982 0.2268 22.483 <2e-16 ***
## poly(Latitude, 3)1 -1.1682 1.1803 -0.990 0.3223
## poly(Latitude, 3)2 -2.3212 1.1633 -1.995 0.0460 *
## poly(Latitude, 3)3 -2.6746 1.1575 -2.311 0.0208 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 -0.067
## ply(Ltt,3)2 0.072 -0.069
## ply(Ltt,3)3 -0.017 0.024 -0.085
r2(Inv_mod3)
## # R2 for Mixed Models
##
## Conditional R2: 0.808
## Marginal R2: 0.437
predictions <- ggpredict(Inv_mod3, terms ="Latitude [all]")
# Plot the GAM model predictions
Inv_plot <- ggplot() +
geom_jitter(data = TransectGroup_wide, aes(x= Latitude, y=Invertebrate_predators), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) +
Theme1 +
labs( y="Invertebrate predators", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y = element_text(size = 7))
Richness_transect$Latitude <- Richness_transect$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
plot(Richness_transect$Latitude, Richness_transect$Richness, xlab = "Latitude", ylab = "Richness", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Richness_transect$Latitude, Richness_transect$Richness, pch = 21, col = mygrey, bg = mygrey)
descdist(Richness_transect$Richness, discrete = FALSE)
## summary statistics
## ------
## min: 2 max: 15
## median: 8
## mean: 8.115385
## estimated sd: 3.076712
## estimated skewness: 0.3378056
## estimated kurtosis: 3.339398
normal_dist <- fitdist(Richness_transect$Richness, "norm")
plot(normal_dist)
R_mod1 <- glmer(Richness ~ Latitude+ (1|Site), data= Richness_transect, family = gaussian())
## Warning in glmer(Richness ~ Latitude + (1 | Site), data = Richness_transect, :
## calling glmer() with family=gaussian (identity link) as a shortcut to lmer() is
## deprecated; please call lmer() directly
R_mod2 <- glmer(Richness ~poly(Latitude,2)+ (1|Site), data= Richness_transect, family = gaussian())
## Warning in glmer(Richness ~ poly(Latitude, 2) + (1 | Site), data =
## Richness_transect, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
R_mod3 <- glmer(Richness ~poly(Latitude,3)+ (1|Site), data= Richness_transect, family = gaussian())
## Warning in glmer(Richness ~ poly(Latitude, 3) + (1 | Site), data =
## Richness_transect, : calling glmer() with family=gaussian (identity link) as a
## shortcut to lmer() is deprecated; please call lmer() directly
AICc(R_mod2, R_mod1, R_mod3)
## df AICc
## R_mod2 5 127.3899
## R_mod1 4 136.8463
## R_mod3 6 124.1031
BIC(R_mod2, R_mod1, R_mod3)
## df BIC
## R_mod2 5 130.6804
## R_mod1 4 139.9739
## R_mod3 6 127.2306
summary(R_mod3)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Richness ~ poly(Latitude, 3) + (1 | Site)
## Data: Richness_transect
##
## REML criterion at convergence: 107.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6252 -0.5226 0.0082 0.3615 2.0845
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.464 2.338
## Residual 4.171 2.042
## Number of obs: 26, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 8.0744 0.7901 10.220
## poly(Latitude, 3)1 -3.4362 4.0846 -0.841
## poly(Latitude, 3)2 -3.0754 4.0715 -0.755
## poly(Latitude, 3)3 -5.9172 4.0749 -1.452
##
## Correlation of Fixed Effects:
## (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 -0.079
## ply(Ltt,3)2 0.057 -0.081
## ply(Ltt,3)3 -0.014 0.039 -0.072
r2(R_mod3)
## # R2 for Mixed Models
##
## Conditional R2: 0.649
## Marginal R2: 0.189
No significant effect
# Plot the GAM model predictions
Richness_plot <- ggplot() +
geom_jitter(data = Richness_transect, aes(x= Latitude, y=Richness), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
Theme1 +
labs( y="Species Richness", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y = element_text(size = 8))
Size_fish <- filter(Size_mean_sp, Type == "Fish")
Size_fish$Latitude <- Size_fish$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
plot(Size_fish$Latitude, Size_fish$mean_Size_sp, xlab = "Latitude", ylab = "Fish Length", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Size_fish$Latitude, Size_fish$mean_Size_sp, pch = 21, col = mygrey, bg = mygrey)
descdist(Size_fish$mean_Size_sp, discrete = FALSE)
## summary statistics
## ------
## min: 5.4 max: 65
## median: 23
## mean: 22.49199
## estimated sd: 9.684538
## estimated skewness: 0.9708081
## estimated kurtosis: 5.717812
normal_dist <- fitdist(Size_fish$mean_Size_sp, "logis")
plot(normal_dist)
Length_mod1 <- glmer(mean_Size_sp ~ Latitude+ (1|Site) + (1|Species), data= Size_fish, family = Gamma(link = "log"))
Length_mod2 <- glmer(mean_Size_sp ~poly(Latitude,2)+ (1|Site) + (1|Species), data= Size_fish, family = Gamma(link = "log"))
Length_mod3 <- glmer(mean_Size_sp ~poly(Latitude,3)+ (1|Site) + (1|Species), data= Size_fish, family = Gamma(link = "log"))
check_model(Length_mod1)
AICc(Length_mod1, Length_mod2, Length_mod3)
## df AICc
## Length_mod1 5 773.3847
## Length_mod2 6 772.3818
## Length_mod3 7 774.6382
BIC(Length_mod1, Length_mod2, Length_mod3)
## df BIC
## Length_mod1 5 786.7493
## Length_mod2 6 788.3065
## Length_mod3 7 793.0830
summary(Length_mod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: mean_Size_sp ~ poly(Latitude, 2) + (1 | Site) + (1 | Species)
## Data: Size_fish
##
## AIC BIC logLik deviance df.resid
## 771.6 788.3 -379.8 759.6 113
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2255 -0.4292 0.0110 0.4239 4.0786
##
## Random effects:
## Groups Name Variance Std.Dev.
## Species (Intercept) 0.091860 0.30308
## Site (Intercept) 0.002929 0.05412
## Residual 0.073020 0.27022
## Number of obs: 119, groups: Species, 28; Site, 12
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 3.0014 0.1174 25.566 <2e-16 ***
## poly(Latitude, 2)1 0.1303 0.3940 0.331 0.7409
## poly(Latitude, 2)2 -0.7157 0.3493 -2.049 0.0405 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(L,2)1
## ply(Ltt,2)1 -0.085
## ply(Ltt,2)2 -0.027 0.036
r2(Length_mod2)
## # R2 for Mixed Models
##
## Conditional R2: 0.585
## Marginal R2: 0.026
predictions <- ggpredict(Length_mod2, terms ="Latitude [all]")
# Plot the GAM model predictions
Size_plot <- ggplot() +
geom_jitter(data = Size_fish, aes(x= Latitude, y=mean_Size_sp), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) +
Theme1 +
labs( y="Fish length (cm)", x="Latitude (°S)") +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(), axis.title.y=element_text(8))
Phlorotannins_metadata$Latitude <- Phlorotannins_metadata$Latitude*-1
# First, examine a simple scatter plot (always a good idea!)
plot(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, xlab = "Latitude", ylab = "phlorotannins", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(Phlorotannins_metadata$Latitude, Phlorotannins_metadata$Phlorotannins, pch = 21, col = mygrey, bg = mygrey)
descdist(Phlorotannins_metadata$Phlorotannins, discrete = FALSE)
## summary statistics
## ------
## min: 1.83 max: 94.69
## median: 11.625
## mean: 17.26067
## estimated sd: 15.2777
## estimated skewness: 2.275762
## estimated kurtosis: 9.890266
normal_dist <- fitdist(Phlorotannins_metadata$Phlorotannins, "exp")
plot(normal_dist)
Phloro_mod1 <- glmer(Phlorotannins ~ Latitude+ (1|Site), data= Phlorotannins_metadata, family = Gamma(link= "log"))
Phloro_mod2 <- glmer(Phlorotannins ~poly(Latitude,2)+ (1|Site), data= Phlorotannins_metadata, family = Gamma(link= "log"))
Phloro_mod3 <- glmer(Phlorotannins ~poly(Latitude,3)+ (1|Site), data= Phlorotannins_metadata, family = Gamma(link= "log"))
check_model(Phloro_mod1)
AICc(Phloro_mod1, Phloro_mod2, Phloro_mod3)
## df AICc
## Phloro_mod1 4 831.2629
## Phloro_mod2 5 832.2302
## Phloro_mod3 6 828.0473
summary(Phloro_mod3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( log )
## Formula: Phlorotannins ~ poly(Latitude, 3) + (1 | Site)
## Data: Phlorotannins_metadata
##
## AIC BIC logLik deviance df.resid
## 827.3 844.0 -407.7 815.3 114
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5999 -0.5788 -0.0462 0.2840 3.8222
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.06263 0.2503
## Residual 0.29291 0.5412
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 2.6600 0.1165 22.838 < 2e-16 ***
## poly(Latitude, 3)1 -3.1963 1.2765 -2.504 0.01228 *
## poly(Latitude, 3)2 -2.0276 1.2753 -1.590 0.11185
## poly(Latitude, 3)3 3.8707 1.2767 3.032 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) p(L,3)1 p(L,3)2
## ply(Ltt,3)1 0.000
## ply(Ltt,3)2 0.001 0.000
## ply(Ltt,3)3 0.000 -0.001 -0.013
r2(Phloro_mod3)
## # R2 for Mixed Models
##
## Conditional R2: 0.546
## Marginal R2: 0.435
predictions <- ggpredict(Phloro_mod3, terms ="Latitude [all]")
# Plot the GAM model predictions
Phloro_plot <- ggplot() +
geom_jitter(data = Phlorotannins_metadata, aes(x= Latitude, y=Phlorotannins), width = 0.08, height = 0.5, size = 1, alpha = 0.6) +
geom_ribbon(data = predictions, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), alpha = 0.2) +
geom_line(data = predictions, aes(x = x, y = predicted), linewidth = 0.6) +
Theme1 +
labs( y="Phlorotannins (mg"~g^-1*" DW)", x="Latitude (°S)")
Temperature_plot + Salinity_plot + Oxygen_plot + Phloro_plot+ plot_layout(ncol=2)
Figure 2: Latitudinal variation of abiotic variables: temperature (°C), salinity (ppt) and dissolved oxygen (mgL) and visibility (m).
Save the plot in tiff format
ggsave("Latitudevsabiotic.jpeg", units="mm", width=125, height=85, dpi=300)
Visibility_plot + Bare_plot+ Kelp_plot + UA_plot + Inv_plot + PredFish_plot + Size_plot + Richness_plot + PCOA_plot+ Herb_plot + plot_layout(ncol=2)
## Warning in grid.Call(C_stringMetric, as.graphicsAnnot(x$label)): font family
## not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y, : font
## family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
## Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, :
## font family not found in Windows font database
Figure 3: Latitudinal variation of biotic and habitat variables
ggsave("FigurelatvsBio.jpeg", units="mm", width=110, height=140, dpi=300)
Latitudinal patterns in consumption (predation and herbivory) were explored and analysed
First change latitude to positive numbers for plotting
Predation_env$Latitude <- Predation_env$Latitude*-1
C_24 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 24)
C_2 <- filter(Predation_env, Bait2 == "Crabs" & Assay_time2 == 2)
S_24 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 24)
S_2 <- filter(Predation_env, Bait2 == "Squidpops" & Assay_time2 == 2)
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(C_2$Latitude, C_2$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(C_2$Latitude, C_2$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)
descdist(C_2$Consumed_individuals, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 5
## median: 1
## mean: 1.958333
## estimated sd: 2.043213
## estimated skewness: 0.4960334
## estimated kurtosis: 1.571122
normal_dist <- fitdist(C_2$Consumed_individuals, "exp")
plot(normal_dist)
C2_mod1 <- glmer(Consumed_individuals ~ Latitude+ (1|Site), data= C_2, family = poisson())
C2_mod2 <- glmer(Consumed_individuals ~poly(Latitude,2)+ (1|Site), data= C_2, family = poisson())
C2_mod3 <- glmer(Consumed_individuals ~poly(Latitude,3)+ (1|Site), data= C_2, family = poisson())
check_model(C2_mod1)
AICc(C2_mod2, C2_mod1, C2_mod3)
## df AICc
## C2_mod2 4 424.9806
## C2_mod1 3 423.1983
## C2_mod3 5 425.7216
BIC(C2_mod2, C2_mod1, C2_mod3)
## df BIC
## C2_mod2 4 435.7828
## C2_mod1 3 431.3539
## C2_mod3 5 439.1328
summary(C2_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Latitude + (1 | Site)
## Data: C_2
##
## AIC BIC logLik deviance df.resid
## 423.0 431.4 -208.5 417.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92403 -0.94180 0.03414 0.62864 3.02889
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.414 0.6434
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.33581 0.91197 2.561 0.0104 *
## Latitude -0.06553 0.03032 -2.161 0.0307 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Latitude -0.975
r2(C2_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.601
## Marginal R2: 0.198
a. Analyses
# First, examine a simple scatter plot (always a good idea!)
plot(C_24$Latitude, C_24$Consumed_individuals, xlab = "Latitude", ylab = "crabs consumed 24 h", las=1)
mygrey <- grey(level = 0.65, alpha = 0.4)
points(C_24$Latitude, C_24$Consumed_individuals, pch = 21, col = mygrey, bg = mygrey)
descdist(C_24$Consumed_individuals, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 5
## median: 5
## mean: 3.791667
## estimated sd: 1.511279
## estimated skewness: -0.8418895
## estimated kurtosis: 2.211133
normal_dist <- fitdist(C_24$Consumed_individuals, "exp")
plot(normal_dist)
C24_mod1 <- glm(Consumed_individuals ~ Latitude, data= C_24, family = poisson())
C24_mod2 <- glm(Consumed_individuals ~poly(Latitude,2), data= C_24, family = poisson())
C24_mod3 <- glm(Consumed_individuals ~poly(Latitude,3), data= C_24, family = poisson())
check_overdispersion(C24_mod1)
## # Overdispersion test
##
## dispersion ratio = 0.420
## Pearson's Chi-Squared = 49.526
## p-value = 1
## No overdispersion detected.
check_model(C24_mod1)
AICc(C24_mod1, C24_mod2, C24_mod3)
## df AICc
## C24_mod1 2 433.7401
## C24_mod2 3 434.2198
## C24_mod3 4 436.0433
BIC(C24_mod1, C24_mod2, C24_mod3)
## df BIC
## C24_mod1 2 439.2125
## C24_mod2 3 442.3754
## C24_mod3 4 446.8454
summary(C24_mod1)
##
## Call:
## glm(formula = Consumed_individuals ~ Latitude, family = poisson(),
## data = C_24)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.452269 0.208192 11.779 < 2e-16 ***
## Latitude -0.038689 0.007222 -5.357 8.46e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 88.137 on 119 degrees of freedom
## Residual deviance: 58.265 on 118 degrees of freedom
## AIC: 433.64
##
## Number of Fisher Scoring iterations: 4
r2(C24_mod1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.424
Calculate predictions for Crab consumption after 2 and 24 h
predictions <- ggpredict(C2_mod1, terms ="Latitude")
predictions1 <- ggpredict(C24_mod1, terms ="Latitude")
Add new column to mix datasets
predictions$Assay_time <- "2 h"
predictions1$Assay_time <- "24 h"
C_2$Assay_time <- "2 h"
C_24$Assay_time <- "24 h"
#Bind datasets
total <- rbind(predictions, predictions1)
total2 <- rbind (C_2, C_24)
ggplot()+
geom_ribbon(data = total, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, group = Assay_time, fill = Assay_time ), alpha = .25) +
geom_line(data = total, aes(x = x, y = predicted, color = Assay_time, group = Assay_time), size = 1) +
geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = 0.5)+
labs(x = "Latitude", y = "Crabs consumed") + Theme1 +
scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
scale_fill_manual("Assay time", values=c("darkslateblue", "orange1")) +
theme(legend.position = c(0.86, 0.81),legend.title = element_text(size = 8),
legend.text = element_text(size = 6.3), legend.background = element_rect(fill="NA",
size=1,
color ="NA")) +
scale_x_continuous(breaks = seq(15, 40, 5)) +
scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 27 rows containing missing values or values outside the scale range
## (`geom_point()`).
Save plot in tiff format
ggsave("Crabs.jpeg", units="mm", width=100, height=65, dpi=300)
## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`geom_point()`).
a. Analyses
descdist(S_2$Consumed_individuals, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 5
## median: 2
## mean: 2.375
## estimated sd: 1.974895
## estimated skewness: 0.1159755
## estimated kurtosis: 1.430103
normal_dist <- fitdist(S_2$Consumed_individuals, "nbinom")
plot(normal_dist)
S2_mod1 <- glmer(Consumed_individuals ~ Latitude + (1|Site), data= S_2, family = poisson())
S2_mod2 <- glmer(Consumed_individuals ~poly(Latitude,2)+ (1|Site), data= S_2, family = poisson())
S2_mod3 <- glmer(Consumed_individuals ~poly(Latitude,3)+ (1|Site), data= S_2, family = poisson())
check_overdispersion(S2_mod1)
## # Overdispersion test
##
## dispersion ratio = 1.050
## Pearson's Chi-Squared = 122.805
## p-value = 0.338
## No overdispersion detected.
check_model(S2_mod1)
AICc(S2_mod2, S2_mod1, S2_mod3)
## df AICc
## S2_mod2 4 462.0740
## S2_mod1 3 460.3905
## S2_mod3 5 463.5723
BIC(S2_mod2, S2_mod1, S2_mod3)
## df BIC
## S2_mod2 4 472.8761
## S2_mod1 3 468.5461
## S2_mod3 5 476.9835
summary(S2_mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Latitude + (1 | Site)
## Data: S_2
##
## AIC BIC logLik deviance df.resid
## 460.2 468.5 -227.1 454.2 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7204 -0.8067 -0.0802 0.7098 2.3659
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5241 0.7239
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.14773 0.99008 0.149 0.881
## Latitude 0.01733 0.03209 0.540 0.589
##
## Correlation of Fixed Effects:
## (Intr)
## Latitude -0.974
r2(S2_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.618
## Marginal R2: 0.016
a. Analyses
descdist(S_24$Consumed_individuals, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 5
## median: 5
## mean: 4.25
## estimated sd: 1.258584
## estimated skewness: -1.822846
## estimated kurtosis: 5.495426
normal_dist <- fitdist(S_24$Consumed_individuals, "pois")
plot(normal_dist)
S24_mod1 <- glm(Consumed_individuals ~ Latitude, data= S_24, family = poisson())
S24_mod2 <- glm(Consumed_individuals ~poly(Latitude,2), data= S_24, family = poisson())
S24_mod3 <- glm(Consumed_individuals ~poly(Latitude,3), data= S_24, family = poisson())
check_overdispersion(S24_mod1)
## # Overdispersion test
##
## dispersion ratio = 0.348
## Pearson's Chi-Squared = 41.044
## p-value = 1
## No overdispersion detected.
check_model(S24_mod1)
AICc(S24_mod2, S24_mod1, S24_mod3)
## df AICc
## S24_mod2 3 447.7243
## S24_mod1 2 446.1367
## S24_mod3 4 449.8652
BIC(S24_mod2, S24_mod1, S24_mod3)
## df BIC
## S24_mod2 3 455.8799
## S24_mod1 2 451.6091
## S24_mod3 4 460.6673
summary(S24_mod1)
##
## Call:
## glm(formula = Consumed_individuals ~ Latitude, family = poisson(),
## data = S_24)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.972899 0.196634 10.033 < 2e-16 ***
## Latitude -0.017890 0.006608 -2.707 0.00678 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 54.280 on 118 degrees of freedom
## AIC: 446.03
##
## Number of Fisher Scoring iterations: 4
r2(S24_mod1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.149
#For squidpops 24 h
predictions3 <- ggpredict(S24_mod1, terms ="Latitude")
Add new column to mix datasets
S_2$Assay_time <- "2 h"
S_24$Assay_time <- "24 h"
predictions3$Assay_time <- "24 h"
##Keep only colums for plotting
total2 <- rbind (S_2, S_24)
ggplot() +
geom_ribbon(data = predictions3, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high, group = Assay_time, fill = Assay_time ), alpha = .25) +
geom_line(data = predictions3, aes(x = x, y = predicted, color = Assay_time, group = Assay_time ), size = 1) +
geom_jitter(data = total2, aes(x= Latitude, y=Consumed_individuals, color =Assay_time), width=0.15, height = 0.2, size = 1, alpha = .5)+
labs(x = "Latitude", y = "Squidpops consumed") + Theme1 +
scale_color_manual("Assay time",values=c("darkslateblue", "orange1"))+
scale_fill_manual("Assay time", values=c("orange1")) +
theme(legend.position = "none") +
scale_x_continuous(breaks = seq(15, 40, 5)) +
scale_y_continuous(limits = c(0,6.5), breaks = seq(0, 5, 1))
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).
Save plot in tiff format
ggsave("Squidpops.jpeg", units="mm", width=100, height=65, dpi=300)
## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
Herbivory_env$Latitude <- Herbivory_env$Latitude*-1
a. Analyses
descdist(Herbivory_env$Consumed_proportion, discrete = FALSE)
## summary statistics
## ------
## min: 0 max: 1
## median: 0
## mean: 0.2260169
## estimated sd: 0.3698909
## estimated skewness: 1.239303
## estimated kurtosis: 2.780775
normal_dist <- fitdist(Herbivory_env$Consumed_proportion, "exp")
plot(normal_dist)
Herb_mod1 <- glmmTMB(Consumed_proportion ~ Latitude+ (1|Site), data= Herbivory_env, family = tweedie())
Herb_mod2 <- glmmTMB(Consumed_proportion ~poly(Latitude,2)+ (1|Site), data= Herbivory_env, family = tweedie())
Herb_mod3 <- glmmTMB(Consumed_proportion ~poly(Latitude,3)+ (1|Site), data= Herbivory_env, family = tweedie())
check_model(Herb_mod1)
## `check_outliers()` does not yet support models of class `glmmTMB`.
AICc(Herb_mod2, Herb_mod1, Herb_mod3)
## df AICc
## Herb_mod2 6 732.3190
## Herb_mod1 5 730.7644
## Herb_mod3 7 732.0831
summary(Herb_mod1)
## Family: tweedie ( log )
## Formula: Consumed_proportion ~ Latitude + (1 | Site)
## Data: Herbivory_env
##
## AIC BIC logLik deviance df.resid
## 730.7 752.6 -360.3 720.7 585
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 6.69 2.586
## Number of obs: 590, groups: Site, 12
##
## Dispersion parameter for tweedie family (): 0.875
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.9565 3.5611 -2.515 0.0119 *
## Latitude 0.1899 0.1138 1.668 0.0953 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(Herb_mod1)
## # R2 for Mixed Models
##
## Conditional R2: 0.956
## Marginal R2: 0.195
b. Plot predicted relationship
ggplot() +
geom_jitter(data = Herbivory_env, aes(x= Latitude, y=Consumed_proportion), color = "orange1", width=0.2, height = 0.02, size = 1, alpha = .5)+
labs(x = "Latitude (°S)", y = "Blade consumed proportion") + Theme1 +
theme(legend.position = "none") +
scale_x_continuous(breaks = seq(15, 40, 5)) +
scale_y_continuous(limits = c(-0.1,1.2), breaks = seq(0, 1, 0.2))
Save plot to tiff format
ggsave("Seaweed24.jpeg", units="mm", width=105, height=65, dpi=300)
Here, the individual effects of each predictor will be tested.
Predation_analysis <- Predation_6 %>%
dplyr::select(-Rope, -Bait, -Camera, -Assay_time, -Initial_Individuals, -Remaining_individuals, -Survival_proportion,-Consumed_proportion)
#Separate by Bait and Assay time
Crabs_24 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 24)
Crabs_2 <- filter(Predation_analysis, Bait2 == "Crabs" & Assay_time2 == 2)
Squidpops_24 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 24)
Squidpops_2 <- filter(Predation_analysis, Bait2 == "Squidpops" & Assay_time2 == 2)
1) Standardise
#Crabs 24h
Crabs24_num <- Crabs_24 %>%
dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude, -Day_length)
Crabs24_stand <- decostand(Crabs24_num, method="range", na.rm = TRUE)
#Test autocorrelation
corrplot(cor(Crabs24_stand), method = "number")
No multicolinearity found
Crabs24_stand$Site <- Crabs_24$Site
Crabs24_stand$Consumed_individuals<- Crabs_24$Consumed_individuals
Crabs24_stand$Latitude<- Crabs_24$Latitude
#Crabs 2h
Crabs2_num <- Crabs_2 %>%
dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude)
Crabs2_stand <- decostand(Crabs2_num, method="range", na.rm = TRUE)
Crabs2_stand$Site <- Crabs_2$Site
Crabs2_stand$Consumed_individuals<- Crabs_2$Consumed_individuals
Crabs2_stand$Latitude<- Crabs_2$Latitude
#Squidpops 24h
Squidpops24_num <- Squidpops_24 %>%
dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude)
Squidpops24_stand <- decostand(Squidpops24_num, method="range", na.rm = TRUE)
Squidpops24_stand$Site <- Squidpops_24$Site
Squidpops24_stand$Consumed_individuals<- Squidpops_24$Consumed_individuals
Squidpops24_stand$Latitude<- Squidpops_24$Latitude
#Squidpops 2h
Squidpops2_num <- Squidpops_2 %>%
dplyr::select(-Site, -Bait2, -Assay_time2, -Consumed_individuals, -Latitude)
Squidpops2_stand <- decostand(Squidpops2_num, method="range", na.rm = TRUE)
Squidpops2_stand$Site <- Squidpops_2$Site
Squidpops2_stand$Consumed_individuals<- Squidpops_2$Consumed_individuals
Squidpops2_stand$Latitude<- Squidpops_2$Latitude
my_glmT1<- glmer(Consumed_individuals~ Temperature + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmT1)
AICc(my_glmT1)
## [1] 420.2918
summary(my_glmT1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 420.1 428.4 -207.0 414.1 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93381 -0.92312 0.06015 0.65176 2.79744
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.2979 0.5458
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8086 0.4487 -1.802 0.07151 .
## Temperature 1.8546 0.6104 3.038 0.00238 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Temperature -0.918
r2(my_glmT1)
## # R2 for Mixed Models
##
## Conditional R2: 0.599
## Marginal R2: 0.307
my_glmSC2_1<- glmer(Consumed_individuals~ Salinity_ppt + (1|Site), data=Crabs2_stand, family = poisson)
check_overdispersion(my_glmSC2_1)
## # Overdispersion test
##
## dispersion ratio = 1.189
## Pearson's Chi-Squared = 139.141
## p-value = 0.08
## No overdispersion detected.
check_model(my_glmSC2_1)
AICc(my_glmSC2_1)
## [1] 424.5386
summary(my_glmSC2_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 424.3 432.7 -209.2 418.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92849 -0.91521 0.04145 0.66421 2.99635
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.468 0.6841
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6810 0.6667 -1.022 0.3070
## Salinity_ppt 1.5821 0.9192 1.721 0.0852 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Salinty_ppt -0.947
r2(my_glmSC2_1)
## # R2 for Mixed Models
##
## Conditional R2: 0.603
## Marginal R2: 0.148
my_glmOxC2_1<- glmer(Consumed_individuals~ D_oxygen_mgL + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmOxC2_1)
AICc(my_glmOxC2_1)
## [1] 422.5227
summary(my_glmOxC2_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 422.3 430.7 -208.2 416.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9665 -0.9020 0.0611 0.5761 3.2095
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.365 0.6042
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.2934 0.3526 -0.832 0.4054
## D_oxygen_mgL 1.5663 0.6451 2.428 0.0152 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## D_oxygn_mgL -0.836
r2(my_glmOxC2_1)
## # R2 for Mixed Models
##
## Conditional R2: 0.585
## Marginal R2: 0.215
my_glmVC2_1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmVC2_1)
AICc(my_glmVC2_1)
## [1] 427.2626
summary(my_glmVC2_1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 427.1 435.4 -210.5 421.1 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93460 -0.89601 0.07717 0.53694 3.02786
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5924 0.7697
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.39200 0.44714 0.877 0.381
## V_Visibility -0.01098 0.89880 -0.012 0.990
##
## Correlation of Fixed Effects:
## (Intr)
## V_Visibilty -0.847
r2(my_glmVC2_1)
## # R2 for Mixed Models
##
## Conditional R2: 0.592
## Marginal R2: 0.000
my_glmC2_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_PCO1)
AICc(my_glmC2_PCO1)
## [1] 426.871
summary(my_glmC2_PCO1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 426.7 435.0 -210.3 420.7 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92572 -0.90531 0.06451 0.54796 2.99751
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5818 0.7628
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5359 0.3320 1.614 0.107
## PCo1 -0.4240 0.6764 -0.627 0.531
##
## Correlation of Fixed Effects:
## (Intr)
## PCo1 -0.703
r2(my_glmC2_PCO1)
## # R2 for Mixed Models
##
## Conditional R2: 0.597
## Marginal R2: 0.023
my_glmC2_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_K1)
AICc(my_glmC2_K1)
## [1] 424.3606
summary(my_glmC2_K1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 424.2 432.5 -209.1 418.2 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92113 -0.89423 0.02569 0.55099 3.16011
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4494 0.6704
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.02724 0.29757 0.092 0.9271
## Kelp 1.32373 0.73207 1.808 0.0706 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Kelp -0.705
r2(my_glmC2_K1)
## # R2 for Mixed Models
##
## Conditional R2: 0.589
## Marginal R2: 0.138
my_glmC2_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_BR1)
AICc(my_glmC2_BR1)
## [1] 424.1326
summary(my_glmC2_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 423.9 432.3 -209.0 417.9 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9360 -0.9040 0.0409 0.5802 3.2149
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4321 0.6573
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8447 0.3094 2.730 0.00633 **
## Bare_rock -1.1994 0.6315 -1.899 0.05751 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Bare_rock -0.741
r2(my_glmC2_BR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.593
## Marginal R2: 0.163
my_glmC2_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_UA1)
AICc(my_glmC2_UA1)
## [1] 427.0992
summary(my_glmC2_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 426.9 435.3 -210.4 420.9 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92743 -0.90039 0.08966 0.54638 3.04101
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5829 0.7635
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4574 0.2908 1.573 0.116
## Understorey_algae -0.2698 0.6637 -0.406 0.684
##
## Correlation of Fixed Effects:
## (Intr)
## Undrstry_lg -0.584
r2(my_glmC2_UA1)
## # R2 for Mixed Models
##
## Conditional R2: 0.592
## Marginal R2: 0.009
my_glmC2_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_SR1)
AICc(my_glmC2_SR1)
## [1] 425.1803
summary(my_glmC2_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 425.0 433.3 -209.5 419.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.90607 -0.90370 0.06207 0.59243 3.08960
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4998 0.707
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1651 0.4393 -0.376 0.707
## Mean_SR 1.1522 0.7741 1.488 0.137
##
## Correlation of Fixed Effects:
## (Intr)
## Mean_SR -0.864
r2(my_glmC2_SR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.598
## Marginal R2: 0.105
my_glmC2_Inv1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_Inv1)
AICc(my_glmC2_Inv1)
## [1] 427.1083
summary(my_glmC2_Inv1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 426.9 435.3 -210.5 420.9 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93478 -0.89816 0.08557 0.54475 3.04770
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5865 0.7658
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2603 0.4017 0.648 0.517
## Invertebrate_predators 0.2879 0.7302 0.394 0.693
##
## Correlation of Fixed Effects:
## (Intr)
## Invrtbrt_pr -0.808
r2(my_glmC2_Inv1)
## # R2 for Mixed Models
##
## Conditional R2: 0.593
## Marginal R2: 0.009
my_glmC2_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_PF1)
AICc(my_glmC2_PF1)
## [1] 426.858
summary(my_glmC2_PF1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 426.7 435.0 -210.3 420.7 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92759 -0.89488 0.07857 0.55013 3.05066
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5723 0.7565
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2205 0.3528 0.625 0.532
## Predatory_fish 0.4397 0.6854 0.642 0.521
##
## Correlation of Fixed Effects:
## (Intr)
## Predtry_fsh -0.748
r2(my_glmC2_PF1)
## # R2 for Mixed Models
##
## Conditional R2: 0.593
## Marginal R2: 0.022
my_glmC2_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Crabs2_stand, family = poisson)
check_model(my_glmC2_FL1)
AICc(my_glmC2_FL1)
## [1] 426.9261
summary(my_glmC2_FL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 426.7 435.1 -210.4 420.7 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93714 -0.91055 0.06116 0.54889 3.03983
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5781 0.7603
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.06379 0.60517 0.105 0.916
## Fish_size 0.52233 0.89451 0.584 0.559
##
## Correlation of Fixed Effects:
## (Intr)
## Fish_size -0.921
r2(my_glmC2_FL1)
## # R2 for Mixed Models
##
## Conditional R2: 0.594
## Marginal R2: 0.020
my_glmC24_T1<- glmer(Consumed_individuals~Temperature + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_T1)
AICc(my_glmC24_T1)
## [1] 438.9369
summary(my_glmC24_T1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 438.7 447.1 -216.4 432.7 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1336 -0.2824 0.1619 0.4568 1.3189
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.0003175 0.01782
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7379 0.1353 5.456 4.88e-08 ***
## Temperature 0.8747 0.1776 4.925 8.42e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Temperature -0.936
r2(my_glmC24_T1)
## # R2 for Mixed Models
##
## Conditional R2: 0.230
## Marginal R2: 0.229
my_glmC24_S1<- glmer(Consumed_individuals~Salinity_ppt + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_S1)
AICc(my_glmC24_S1)
## [1] 443.1697
summary(my_glmC24_S1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 443.0 451.3 -218.5 437.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.01001 -0.46529 0.09323 0.42898 1.47247
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.01181 0.1087
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.6358 0.1945 3.269 0.001078 **
## Salinity_ppt 0.9897 0.2621 3.777 0.000159 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Salinty_ppt -0.956
r2(my_glmC24_S1)
## # R2 for Mixed Models
##
## Conditional R2: 0.233
## Marginal R2: 0.195
my_glmC24_Ox1<- glmer(Consumed_individuals~D_oxygen_mgL + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_Ox1)
AICc(my_glmC24_Ox1)
## [1] 452.5543
summary(my_glmC24_Ox1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 452.3 460.7 -223.2 446.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0829 -0.3910 0.1934 0.3175 1.2599
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.05924 0.2434
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1962 0.1534 7.796 6.38e-15 ***
## D_oxygen_mgL 0.2413 0.2866 0.842 0.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## D_oxygn_mgL -0.831
r2(my_glmC24_Ox1)
## # R2 for Mixed Models
##
## Conditional R2: 0.216
## Marginal R2: 0.017
my_glmC24_V1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_V1)
AICc(my_glmC24_V1)
## [1] 452.5114
summary(my_glmC24_V1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 452.3 460.7 -223.2 446.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0495 -0.3868 0.2021 0.3902 1.2724
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.05948 0.2439
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4210 0.1604 8.858 <2e-16 ***
## V_Visibility -0.2825 0.3267 -0.865 0.387
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## V_Visibilty -0.846
r2(my_glmC24_V1)
## # R2 for Mixed Models
##
## Conditional R2: 0.218
## Marginal R2: 0.019
my_glmC24_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_PCO1)
AICc(my_glmC24_PCO1)
## [1] 450.5558
summary(my_glmC24_PCO1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 450.3 458.7 -222.2 444.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9758 -0.3167 0.2087 0.3774 1.2274
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.04741 0.2177
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4411 0.1109 12.993 <2e-16 ***
## PCo1 -0.3920 0.2300 -1.705 0.0883 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PCo1 -0.698
r2(my_glmC24_PCO1)
## # R2 for Mixed Models
##
## Conditional R2: 0.223
## Marginal R2: 0.065
my_glmC24_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_K1)
AICc(my_glmC24_K1)
## [1] 448.7633
summary(my_glmC24_K1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 448.6 456.9 -221.3 442.6 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0196 -0.4200 0.1414 0.4214 1.2459
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.03546 0.1883
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1427 0.1040 10.987 <2e-16 ***
## Kelp 0.5844 0.2518 2.321 0.0203 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Kelp -0.715
r2(my_glmC24_K1)
## # R2 for Mixed Models
##
## Conditional R2: 0.210
## Marginal R2: 0.090
my_glmC24_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_BR1)
AICc(my_glmC24_BR1)
## [1] 447.5098
summary(my_glmC24_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 447.3 455.7 -220.7 441.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0745 -0.3637 0.1577 0.3357 1.5700
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.02761 0.1662
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.51816 0.09993 15.192 < 2e-16 ***
## Bare_rock -0.56632 0.20938 -2.705 0.00683 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Bare_rock -0.733
r2(my_glmC24_BR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.216
## Marginal R2: 0.123
my_glmC24_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_UA1)
AICc(my_glmC24_UA1)
## [1] 453.1994
summary(my_glmC24_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 453.0 461.4 -223.5 447.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0561 -0.4182 0.1760 0.3756 1.3128
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.06408 0.2531
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.28831 0.10876 11.846 <2e-16 ***
## Understorey_algae 0.05185 0.24524 0.211 0.833
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Undrstry_lg -0.590
r2(my_glmC24_UA1)
## # R2 for Mixed Models
##
## Conditional R2: 0.216
## Marginal R2: 0.001
my_glmC24_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_SR1)
AICc(my_glmC24_SR1)
## [1] 451.6355
summary(my_glmC24_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 451.4 459.8 -222.7 445.4 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9950 -0.3586 0.1690 0.4479 1.2775
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.05277 0.2297
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1192 0.1647 6.797 1.07e-11 ***
## Mean_SR 0.3830 0.2925 1.309 0.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mean_SR -0.867
r2(my_glmC24_SR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.216
## Marginal R2: 0.040
my_glmC24_Inv1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_Inv1)
AICc(my_glmC24_Inv1)
## [1] 453.0005
summary(my_glmC24_Inv1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 452.8 461.2 -223.4 446.8 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0448 -0.4524 0.1776 0.4005 1.3483
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.06186 0.2487
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3610 0.1462 9.308 <2e-16 ***
## Invertebrate_predators -0.1345 0.2703 -0.498 0.619
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Invrtbrt_pr -0.805
r2(my_glmC24_Inv1)
## # R2 for Mixed Models
##
## Conditional R2: 0.214
## Marginal R2: 0.007
my_glmC24_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_PF1)
AICc(my_glmC24_PF1)
## [1] 452.9183
summary(my_glmC24_PF1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 452.7 461.1 -223.4 446.7 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0344 -0.3964 0.1898 0.4233 1.2267
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.06197 0.2489
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2463 0.1308 9.528 <2e-16 ***
## Predatory_fish 0.1465 0.2549 0.575 0.565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Predtry_fsh -0.748
r2(my_glmC24_PF1)
## # R2 for Mixed Models
##
## Conditional R2: 0.216
## Marginal R2: 0.008
my_glmC24_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Crabs24_stand, family = poisson)
check_model(my_glmC24_FL1)
AICc(my_glmC24_FL1)
## [1] 448.4084
summary(my_glmC24_FL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 448.2 456.6 -221.1 442.2 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0671 -0.2882 0.2492 0.3076 1.0942
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.03566 0.1888
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.8921 0.1932 4.618 3.88e-06 ***
## Fish_size 0.6622 0.2796 2.369 0.0179 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Fish_size -0.926
r2(my_glmC24_FL1)
## # R2 for Mixed Models
##
## Conditional R2: 0.226
## Marginal R2: 0.108
my_glmS2_T1<- glmer(Consumed_individuals~Temperature + (1|Site), data=Squidpops2_stand, family = poisson)
check_overdispersion(my_glmS2_T1)
## # Overdispersion test
##
## dispersion ratio = 1.054
## Pearson's Chi-Squared = 123.348
## p-value = 0.326
## No overdispersion detected.
check_model(my_glmS2_T1)
AICc(my_glmS2_T1)
## [1] 460.4034
summary(my_glmS2_T1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Temperature + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.2 468.6 -227.1 454.2 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7049 -0.8189 -0.1183 0.6789 2.3417
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5106 0.7145
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4209 0.5139 0.819 0.413
## Temperature 0.3832 0.7173 0.534 0.593
##
## Correlation of Fixed Effects:
## (Intr)
## Temperature -0.903
r2(my_glmS2_T1)
## # R2 for Mixed Models
##
## Conditional R2: 0.612
## Marginal R2: 0.016
my_glmS2_S1<- glmer(Consumed_individuals~Salinity_ppt + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_S1)
AICc(my_glmS2_S1)
## [1] 460.4596
summary(my_glmS2_S1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Salinity_ppt + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.3 468.6 -227.1 454.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72415 -0.80815 -0.08599 0.69023 2.37161
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5202 0.7212
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9474 0.6296 1.505 0.132
## Salinity_ppt -0.4192 0.8834 -0.475 0.635
##
## Correlation of Fixed Effects:
## (Intr)
## Salinty_ppt -0.935
r2(my_glmS2_S1)
## # R2 for Mixed Models
##
## Conditional R2: 0.615
## Marginal R2: 0.012
my_glmS2_Ox1<- glmer(Consumed_individuals~D_oxygen_mgL + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_Ox1)
AICc(my_glmS2_Ox1)
## [1] 457.1796
summary(my_glmS2_Ox1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ D_oxygen_mgL + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 457.0 465.3 -225.5 451.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72360 -0.79117 -0.07451 0.73700 2.40980
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.3718 0.6098
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1137 0.3581 0.318 0.7508
## D_oxygen_mgL 1.2713 0.6529 1.947 0.0515 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## D_oxygn_mgL -0.844
r2(my_glmS2_Ox1)
## # R2 for Mixed Models
##
## Conditional R2: 0.606
## Marginal R2: 0.166
my_glmS2_V1<- glmer(Consumed_individuals~V_Visibility + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_V1)
AICc(my_glmS2_V1)
## [1] 460.1753
summary(my_glmS2_V1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ V_Visibility + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.0 468.3 -227.0 454.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7021 -0.8107 -0.1061 0.7043 2.4274
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5109 0.7148
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4194 0.4157 1.009 0.313
## V_Visibility 0.5814 0.8141 0.714 0.475
##
## Correlation of Fixed Effects:
## (Intr)
## V_Visibilty -0.847
r2(my_glmS2_V1)
## # R2 for Mixed Models
##
## Conditional R2: 0.617
## Marginal R2: 0.028
my_glmS2_PCO1<- glmer(Consumed_individuals~PCo1 + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_PCO1)
AICc(my_glmS2_PCO1)
## [1] 457.9922
summary(my_glmS2_PCO1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ PCo1 + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 457.8 466.1 -225.9 451.8 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.74367 -0.78077 -0.08552 0.71274 2.46497
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.3944 0.628
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3370 0.2866 1.176 0.2396
## PCo1 0.9343 0.5457 1.712 0.0869 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PCo1 -0.726
r2(my_glmS2_PCO1)
## # R2 for Mixed Models
##
## Conditional R2: 0.604
## Marginal R2: 0.133
my_glmS2_K1<- glmer(Consumed_individuals~Kelp + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_K1)
AICc(my_glmS2_K1)
## [1] 459.4048
summary(my_glmS2_K1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Kelp + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 459.2 467.6 -226.6 453.2 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6926 -0.8389 -0.1096 0.6556 2.4268
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4541 0.6739
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.4384 0.2948 1.487 0.137
## Kelp 0.8480 0.7294 1.163 0.245
##
## Correlation of Fixed Effects:
## (Intr)
## Kelp -0.702
r2(my_glmS2_K1)
## # R2 for Mixed Models
##
## Conditional R2: 0.606
## Marginal R2: 0.067
my_glmS2_BR1<- glmer(Consumed_individuals~Bare_rock + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_BR1)
AICc(my_glmS2_BR1)
## [1] 460.4304
summary(my_glmS2_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Bare_rock + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.2 468.6 -227.1 454.2 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7003 -0.8230 -0.1211 0.6796 2.3980
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5128 0.7161
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7914 0.3284 2.410 0.016 *
## Bare_rock -0.3255 0.6415 -0.507 0.612
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Bare_rock -0.738
r2(my_glmS2_BR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.612
## Marginal R2: 0.014
my_glmS2_UA1<- glmer(Consumed_individuals~Understorey_algae + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_UA1)
AICc(my_glmS2_UA1)
## [1] 460.6711
summary(my_glmS2_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Understorey_algae + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.5 468.8 -227.2 454.5 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.71045 -0.81805 -0.09991 0.68648 2.36196
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5282 0.7267
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.68441 0.27845 2.458 0.014 *
## Understorey_algae -0.06958 0.62237 -0.112 0.911
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Undrstry_lg -0.592
r2(my_glmS2_UA1)
## # R2 for Mixed Models
##
## Conditional R2: 0.614
## Marginal R2: 0.001
my_glmS2_SR1<- glmer(Consumed_individuals~Mean_SR + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_SR1)
AICc(my_glmS2_SR1)
## [1] 460.6386
summary(my_glmS2_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Mean_SR + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.4 468.8 -227.2 454.4 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.71377 -0.81473 -0.09892 0.68827 2.36015
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5293 0.7275
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.7440 0.4298 1.731 0.0835 .
## Mean_SR -0.1641 0.7736 -0.212 0.8320
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mean_SR -0.853
r2(my_glmS2_SR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.615
## Marginal R2: 0.003
my_glmS2_I1<- glmer(Consumed_individuals~Invertebrate_predators + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_I1)
AICc(my_glmS2_I1)
## [1] 460.0393
summary(my_glmS2_I1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Invertebrate_predators + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 459.8 468.2 -226.9 453.8 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7185 -0.8250 -0.1281 0.7024 2.4144
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4876 0.6983
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.9020 0.3548 2.542 0.011 *
## Invertebrate_predators -0.5317 0.6507 -0.817 0.414
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Invrtbrt_pr -0.792
r2(my_glmS2_I1)
## # R2 for Mixed Models
##
## Conditional R2: 0.609
## Marginal R2: 0.036
my_glmS2_PF1<- glmer(Consumed_individuals~Predatory_fish + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_PF1)
AICc(my_glmS2_PF1)
## [1] 460.5443
summary(my_glmS2_PF1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Predatory_fish + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 460.3 468.7 -227.2 454.3 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.70479 -0.82195 -0.09594 0.67244 2.34183
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.5245 0.7242
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.5755 0.3320 1.733 0.083 .
## Predatory_fish 0.2390 0.6392 0.374 0.709
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Predtry_fsh -0.739
r2(my_glmS2_PF1)
## # R2 for Mixed Models
##
## Conditional R2: 0.615
## Marginal R2: 0.008
my_glmS2_FL1<- glmer(Consumed_individuals~Fish_size + (1|Site), data=Squidpops2_stand, family = poisson)
check_model(my_glmS2_FL1)
AICc(my_glmS2_FL1)
## [1] 458.2359
summary(my_glmS2_FL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Fish_size + (1 | Site)
## Data: Squidpops2_stand
##
## AIC BIC logLik deviance df.resid
## 458.0 466.4 -226.0 452.0 117
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72996 -0.80127 -0.09444 0.71113 2.33612
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.4281 0.6543
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.4069 0.4953 2.841 0.0045 **
## Fish_size -1.2026 0.7565 -1.590 0.1119
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Fish_size -0.911
r2(my_glmS2_FL1)
## # R2 for Mixed Models
##
## Conditional R2: 0.617
## Marginal R2: 0.124
my_glmS24_T1<- glmer(Consumed_individuals~Temperature+ (1|Site), data=Squidpops24_stand, family = poisson)
## boundary (singular) fit: see help('isSingular')
my_glmS24_T1_v2<- glm(Consumed_individuals~Temperature, data=Squidpops24_stand, family = poisson)
check_overdispersion(my_glmS24_T1_v2)
## # Overdispersion test
##
## dispersion ratio = 0.331
## Pearson's Chi-Squared = 39.057
## p-value = 1
## No overdispersion detected.
check_model(my_glmS24_T1_v2)
AICc(my_glmS24_T1_v2)
## [1] 443.5292
summary(my_glmS24_T1_v2)
##
## Call:
## glm(formula = Consumed_individuals ~ Temperature, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1227 0.1166 9.632 < 2e-16 ***
## Temperature 0.4879 0.1575 3.098 0.00195 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 51.673 on 118 degrees of freedom
## AIC: 443.43
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_T1_v2)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.200
my_glmS24_S1<- glm(Consumed_individuals~Salinity_ppt, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_S1)
AICc(my_glmS24_S1)
## [1] 447.5637
summary(my_glmS24_S1)
##
## Call:
## glm(formula = Consumed_individuals ~ Salinity_ppt, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1271 0.1434 7.857 3.93e-15 ***
## Salinity_ppt 0.4670 0.1956 2.387 0.017 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 55.707 on 118 degrees of freedom
## AIC: 447.46
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_S1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.122
my_glmS24_O1<- glm(Consumed_individuals~D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_O1)
AICc(my_glmS24_O1)
## [1] 452.8848
summary(my_glmS24_O1)
##
## Call:
## glm(formula = Consumed_individuals ~ D_oxygen_mgL, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.39170 0.08033 17.324 <2e-16 ***
## D_oxygen_mgL 0.12467 0.14949 0.834 0.404
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 61.028 on 118 degrees of freedom
## AIC: 452.78
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_O1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.014
my_glmS24_V1<- glm(Consumed_individuals~V_Visibility, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_V1)
AICc(my_glmS24_V1)
## [1] 452.1049
summary(my_glmS24_V1)
##
## Call:
## glm(formula = Consumed_individuals ~ V_Visibility, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.53164 0.08209 18.659 <2e-16 ***
## V_Visibility -0.20381 0.16916 -1.205 0.228
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 60.248 on 118 degrees of freedom
## AIC: 452
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_V1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.030
my_glmS24_PC1<- glm(Consumed_individuals~PCo1, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_PC1)
AICc(my_glmS24_PC1)
## [1] 451.8008
summary(my_glmS24_PC1)
##
## Call:
## glm(formula = Consumed_individuals ~ PCo1, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.50543 0.06171 24.394 <2e-16 ***
## PCo1 -0.16872 0.12776 -1.321 0.187
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 59.944 on 118 degrees of freedom
## AIC: 451.7
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_PC1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.037
my_glmS24_K1<- glm(Consumed_individuals~Kelp, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_K1)
AICc(my_glmS24_K1)
## [1] 450.7603
summary(my_glmS24_K1)
##
## Call:
## glm(formula = Consumed_individuals ~ Kelp, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.37365 0.06285 21.856 <2e-16 ***
## Kelp 0.25829 0.15153 1.705 0.0883 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 58.904 on 118 degrees of freedom
## AIC: 450.66
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_K1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.058
my_glmS24_BR1<- glm(Consumed_individuals~Bare_rock, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_BR1)
AICc(my_glmS24_BR1)
## [1] 451.3986
summary(my_glmS24_BR1)
##
## Call:
## glm(formula = Consumed_individuals ~ Bare_rock, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.51985 0.06556 23.183 <2e-16 ***
## Bare_rock -0.19684 0.13441 -1.464 0.143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 59.542 on 118 degrees of freedom
## AIC: 451.3
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_BR1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.045
my_glmS24_UA1<- glm(Consumed_individuals~Understorey_algae, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_UA1)
AICc(my_glmS24_UA1)
## [1] 453.5612
summary(my_glmS24_UA1)
##
## Call:
## glm(formula = Consumed_individuals ~ Understorey_algae, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.44308 0.05495 26.261 <2e-16 ***
## Understorey_algae 0.01475 0.12466 0.118 0.906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 61.705 on 118 degrees of freedom
## AIC: 453.46
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_UA1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.000
my_glmS24_SR1<- glm(Consumed_individuals~Mean_SR, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_SR1)
AICc(my_glmS24_SR1)
## [1] 453.2822
summary( my_glmS24_SR1)
##
## Call:
## glm(formula = Consumed_individuals ~ Mean_SR, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.40627 0.08759 16.054 <2e-16 ***
## Mean_SR 0.08470 0.15635 0.542 0.588
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 61.426 on 118 degrees of freedom
## AIC: 453.18
##
## Number of Fisher Scoring iterations: 4
r2( my_glmS24_SR1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.006
my_glmS24_I1<- glm(Consumed_individuals~Invertebrate_predators, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_I1)
AICc(my_glmS24_I1)
## [1] 453.4745
summary(my_glmS24_I1)
##
## Call:
## glm(formula = Consumed_individuals ~ Invertebrate_predators,
## family = poisson, data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.46575 0.07385 19.848 <2e-16 ***
## Invertebrate_predators -0.04313 0.13606 -0.317 0.751
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 61.618 on 118 degrees of freedom
## AIC: 453.37
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_I1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.002
my_glmS24_PF1<- glm(Consumed_individuals~Predatory_fish, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_PF1)
AICc(my_glmS24_PF1)
## [1] 451.6715
summary(my_glmS24_PF1)
##
## Call:
## glm(formula = Consumed_individuals ~ Predatory_fish, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.37803 0.06749 20.418 <2e-16 ***
## Predatory_fish 0.17727 0.12762 1.389 0.165
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 59.815 on 118 degrees of freedom
## AIC: 451.57
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_PF1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.039
my_glmS24_FL1<- glm(Consumed_individuals~Fish_size, data=Squidpops24_stand, family = poisson)
check_model(my_glmS24_FL1)
AICc(my_glmS24_FL1)
## [1] 451.6657
summary(my_glmS24_FL1)
##
## Call:
## glm(formula = Consumed_individuals ~ Fish_size, family = poisson,
## data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.3041 0.1144 11.400 <2e-16 ***
## Fish_size 0.2283 0.1665 1.371 0.17
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 59.809 on 118 degrees of freedom
## AIC: 451.56
##
## Number of Fisher Scoring iterations: 4
r2(my_glmS24_FL1)
## # R2 for Generalized Linear Regression
## Nagelkerke's R2: 0.039
my_glm<- glm(Consumed_individuals~ Temperature * D_oxygen_mgL, data=Squidpops24_stand, family = poisson)
summary(my_glm)
##
## Call:
## glm(formula = Consumed_individuals ~ Temperature * D_oxygen_mgL,
## family = poisson, data = Squidpops24_stand)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.2170 0.1628 7.474 7.78e-14 ***
## Temperature 0.4416 0.2252 1.961 0.0499 *
## D_oxygen_mgL -0.4196 0.4782 -0.877 0.3803
## Temperature:D_oxygen_mgL 0.3665 0.5390 0.680 0.4965
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 61.719 on 119 degrees of freedom
## Residual deviance: 50.734 on 116 degrees of freedom
## AIC: 446.49
##
## Number of Fisher Scoring iterations: 4
AICc(my_glm)
## [1] 446.8355
First, create datasets to be used in the analysis
Herbivory_analysis <- Herbivory_6 %>%
dplyr::select(-Rope, -Bait, -Camera, -Hours, -Initial_length_mm, -Final_Length_mm, -Scraped, -Survival_proportion)
1) Standardise
Herb_all_num <- Herbivory_analysis %>%
dplyr::select(-Site, -Bait2,-Consumed_proportion, -Day_length)
Herb_all_stand <- decostand(Herb_all_num, method="range")
#Test autocorrelation
corrplot(cor(Herb_all_stand), method = "number")
Herb_all_stand$Site <- Herbivory_analysis$Site
Herb_all_stand$Consumption<- Herbivory_analysis$Consumption
Herb_all_stand$Consumed_proportion<- Herbivory_analysis$Consumed_proportion
my_glmL_T1<- glmer(Consumed_proportion~ Temperature + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_overdispersion(my_glmL_T1)
## # Overdispersion test
##
## dispersion ratio = 0.876
## p-value = 0.584
## No overdispersion detected.
check_model(my_glmL_T1)
AICc(my_glmL_T1)
## [1] 444.0932
summary(my_glmL_T1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Temperature + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 444.1 457.2 -219.0 438.1 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87459 -0.28478 -0.08928 0.12354 3.01044
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 6.183 2.487
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.962 1.740 -0.553 0.58
## Temperature -2.399 2.515 -0.954 0.34
##
## Correlation of Fixed Effects:
## (Intr)
## Temperature -0.883
r2(my_glmL_T1)
## # R2 for Mixed Models
##
## Conditional R2: 0.671
## Marginal R2: 0.052
my_glmL_S1<- glmer(Consumed_proportion~ Salinity_ppt + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_S1)
AICc(my_glmL_S1)
## [1] 441.0416
summary(my_glmL_S1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Salinity_ppt + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 441.0 454.1 -217.5 435.0 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.93354 -0.27675 -0.09387 0.13458 2.96049
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 4.546 2.132
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.288 1.846 0.698 0.4853
## Salinity_ppt -5.561 2.705 -2.055 0.0398 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Salinty_ppt -0.922
r2(my_glmL_S1)
## # R2 for Mixed Models
##
## Conditional R2: 0.662
## Marginal R2: 0.194
my_glmL_O1<- glmer(Consumed_proportion~ D_oxygen_mgL + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_O1)
AICc(my_glmL_O1)
## [1] 444.2172
summary(my_glmL_O1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ D_oxygen_mgL + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 444.2 457.3 -219.1 438.2 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.86805 -0.26809 -0.09095 0.11924 2.97058
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 6.02 2.454
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.520 1.451 -2.425 0.0153 *
## D_oxygen_mgL 2.400 2.654 0.904 0.3659
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## D_oxygn_mgL -0.831
r2(my_glmL_O1)
## # R2 for Mixed Models
##
## Conditional R2: 0.664
## Marginal R2: 0.050
my_glmL_V1<- glmer(Consumed_proportion~ V_Visibility + (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_V1)
AICc(my_glmL_V1)
## [1] 444.6448
summary(my_glmL_V1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ V_Visibility + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 444.6 457.7 -219.3 438.6 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.86293 -0.28102 -0.08826 0.12687 3.02234
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 6.441 2.538
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.223 1.509 -2.136 0.0327 *
## V_Visibility 1.713 2.851 0.601 0.5480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## V_Visibilty -0.835
r2(my_glmL_V1)
## # R2 for Mixed Models
##
## Conditional R2: 0.669
## Marginal R2: 0.021
my_glmL_PC1<- glmer(Consumed_proportion~ PCo1+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_PC1)
AICc(my_glmL_PC1)
## [1] 444.2705
summary(my_glmL_PC1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ PCo1 + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 444.2 457.4 -219.1 438.2 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88774 -0.28227 -0.08271 0.12286 3.01571
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 6.24 2.498
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.159 1.175 -2.687 0.00721 **
## PCo1 1.863 2.169 0.859 0.39060
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## PCo1 -0.715
r2(my_glmL_PC1)
## # R2 for Mixed Models
##
## Conditional R2: 0.670
## Marginal R2: 0.044
my_glmL_K1<- glmmTMB(Consumed_proportion~Kelp+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
check_model(my_glmL_K1)
## `check_outliers()` does not yet support models of class `glmmTMB`.
AICc(my_glmL_K1)
## [1] 451.4504
Linear model failed to converge, quadratic model was chosen for comparisons.
summary(my_glmL_K1)
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Kelp + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 451.4 464.5 -222.7 445.4 587
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.716 2.391
## Number of obs: 590, groups: Site, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1796 1.0412 -2.093 0.0363 *
## Kelp -0.8131 2.5535 -0.318 0.7502
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(my_glmL_K1)
## # R2 for Mixed Models
##
## Conditional R2: 0.637
## Marginal R2: 0.006
my_glmL_BR1<- glmer(Consumed_proportion~Bare_rock+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_BR1)
AICc(my_glmL_BR1)
## [1] 443.0529
summary(my_glmL_BR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Bare_rock + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 443.0 456.2 -218.5 437.0 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.88882 -0.26670 -0.09483 0.12657 2.98644
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.49 2.343
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.630 1.189 -3.054 0.00226 **
## Bare_rock 3.035 2.136 1.421 0.15545
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Bare_rock -0.758
r2(my_glmL_BR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.665
## Marginal R2: 0.107
my_glmL_UA1<- glmer(Consumed_proportion~Understorey_algae+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_UA1)
AICc(my_glmL_UA1)
## [1] 437.193
summary(my_glmL_UA1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Understorey_algae + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 437.2 450.3 -215.6 431.2 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8504 -0.2511 -0.0863 0.1267 5.0475
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 3.436 1.854
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0496 0.7058 -1.487 0.13698
## Understorey_algae -6.1317 2.3255 -2.637 0.00837 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Undrstry_lg -0.445
r2(my_glmL_UA1)
## # R2 for Mixed Models
##
## Conditional R2: 0.714
## Marginal R2: 0.415
my_glmL_SR1<- glmer(Consumed_proportion~Mean_SR+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_SR1)
AICc(my_glmL_SR1)
## [1] 444.5391
summary(my_glmL_SR1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Mean_SR + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 444.5 457.6 -219.2 438.5 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.85322 -0.27247 -0.09101 0.12793 2.97373
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 6.295 2.509
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.364 1.556 -2.162 0.0306 *
## Mean_SR 1.834 2.674 0.686 0.4926
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Mean_SR -0.849
r2(my_glmL_SR1)
## # R2 for Mixed Models
##
## Conditional R2: 0.666
## Marginal R2: 0.027
my_glmL_H1<- glmmTMB(Consumed_proportion~Herbivores+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
check_model(my_glmL_H1)
## `check_outliers()` does not yet support models of class `glmmTMB`.
AICc(my_glmL_H1)
## [1] 451.3434
summary(my_glmL_H1)
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ Herbivores + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 451.3 464.4 -222.7 445.3 587
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## Site (Intercept) 5.732 2.394
## Number of obs: 590, groups: Site, 12
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6513 0.9626 -2.754 0.00588 **
## Herbivores 1.1162 2.4550 0.455 0.64934
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2(my_glmL_H1)
## # R2 for Mixed Models
##
## Conditional R2: 0.640
## Marginal R2: 0.012
my_glmL_PL1<- glmer(Consumed_proportion~phloro_mean+ (1|Site), data=Herb_all_stand, family = binomial())
## Warning in eval(family$initialize, rho): non-integer #successes in a binomial
## glm!
check_model(my_glmL_PL1)
AICc(my_glmL_PL1)
## [1] 431.1052
summary(my_glmL_PL1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: Consumed_proportion ~ phloro_mean + (1 | Site)
## Data: Herb_all_stand
##
## AIC BIC logLik deviance df.resid
## 431.1 444.2 -212.5 425.1 587
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8432 -0.3082 -0.0483 0.1339 27.0033
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 2.507 1.583
## Number of obs: 590, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.08397 0.76929 0.109 0.91309
## phloro_mean -9.85455 3.39879 -2.899 0.00374 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## phloro_mean -0.689
r2(my_glmL_PL1)
## # R2 for Mixed Models
##
## Conditional R2: 0.792
## Marginal R2: 0.633
my_glm<- glmer(Consumed_individuals~ Temperature*D_oxygen_mgL + (1|Site), data=Crabs2_stand, family = poisson())
summary(my_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Temperature * D_oxygen_mgL + (1 | Site)
## Data: Crabs2_stand
##
## AIC BIC logLik deviance df.resid
## 421.4 435.3 -205.7 411.4 115
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96584 -0.90387 0.06297 0.57760 2.96316
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.2205 0.4696
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.6119 0.6037 -1.014 0.311
## Temperature 1.0211 0.8251 1.238 0.216
## D_oxygen_mgL -0.2223 1.7120 -0.130 0.897
## Temperature:D_oxygen_mgL 1.3775 1.9254 0.715 0.474
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr D_xy_L
## Temperature -0.849
## D_oxygn_mgL -0.735 0.519
## Tmprtr:D__L 0.730 -0.674 -0.940
AICc(my_glm)
## [1] 421.9314
my_glm<- glmer(Consumed_individuals~ Temperature*Salinity_ppt + (1|Site), data=Crabs24_stand, family = poisson())
## boundary (singular) fit: see help('isSingular')
summary(my_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Temperature * Salinity_ppt + (1 | Site)
## Data: Crabs24_stand
##
## AIC BIC logLik deviance df.resid
## 435.8 449.7 -212.9 425.8 115
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1075 -0.2446 0.1351 0.2772 1.3584
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0 0
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.3196 0.6541 0.489 0.625
## Temperature 0.8991 1.2032 0.747 0.455
## Salinity_ppt 0.8443 1.0544 0.801 0.423
## Temperature:Salinity_ppt -0.3734 1.7801 -0.210 0.834
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr Slnty_
## Temperature -0.964
## Salinty_ppt -0.978 0.944
## Tmprtr:Sln_ 0.961 -0.987 -0.972
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
AICc(my_glm)
## [1] 436.2851
my_glm<- glmer(Consumed_individuals~ Temperature*Salinity_ppt + (1|Site), data=Squidpops24_stand, family = poisson())
## boundary (singular) fit: see help('isSingular')
summary(my_glm)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Consumed_individuals ~ Temperature * Salinity_ppt + (1 | Site)
## Data: Squidpops24_stand
##
## AIC BIC logLik deviance df.resid
## 448 462 -219 438 115
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8095 -0.1049 0.1330 0.1700 0.9803
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0 0
## Number of obs: 120, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.99054 0.58952 1.680 0.0929 .
## Temperature 0.43957 1.13231 0.388 0.6979
## Salinity_ppt 0.28427 0.95104 0.299 0.7650
## Temperature:Salinity_ppt -0.06261 1.66684 -0.038 0.9700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr Slnty_
## Temperature -0.969
## Salinty_ppt -0.979 0.948
## Tmprtr:Sln_ 0.966 -0.988 -0.974
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
AICc(my_glm)
## [1] 448.5681
Make predictions according to individual models with lowest AICc
p_C2T <- ggpredict(my_glmT1, terms ="Temperature [all]")
p_C2Ox <- ggpredict(my_glmOxC2_1, terms ="D_oxygen_mgL [all]")
p_C24T <- ggpredict(my_glmC24_T1, terms ="Temperature [all]")
p_C24S <- ggpredict(my_glmC24_S1, terms ="Salinity_ppt [all]")
p_C24BR <- ggpredict(my_glmC24_BR1, terms ="Bare_rock [all]")
p_C24K <- ggpredict(my_glmC24_K1, terms ="Kelp [all]")
p_C24FL <- ggpredict(my_glmC24_FL1, terms ="Fish_size [all]")
1) Temperature vs Crabs
ggplot() +
geom_ribbon(data = p_C2T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "darkslateblue", alpha = .25) +
geom_line(data = p_C2T, aes(x = x, y =predicted), color = "darkslateblue", linewidth = 0.6) +
geom_jitter(data = Crabs2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.008, height = 0.1, size = 1, alpha = 0.5)+
geom_ribbon(data = p_C24T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_C24T, aes(x = x, y =predicted), color = "orange1", linewidth = 0.6) +
geom_jitter(data = Crabs24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.008, height = 0.1, size = 1, alpha = 0.5)+
labs(x = "Temperature (°C)", y = "Crabs consumed") + Theme1 +
scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1))
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("TempC.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
2) Salinity vs Crabs
ggplot() +
geom_ribbon(data = p_C24S, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_C24S, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Crabs24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
geom_jitter(data = Crabs2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
labs(x = "Salinity (ppt)", y = "Crabs consumed") + Theme1 +
scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,8), breaks = seq(0, 6, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("SalC.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
3) Oxygen vs Crabs
ggplot() +
geom_ribbon(data = p_C2Ox, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "darkslateblue", alpha = .25) +
geom_line(data = p_C2Ox, aes(x = x, y = predicted), color = "darkslateblue", size = 1) +
geom_jitter(data = Crabs2_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
geom_jitter(data = Crabs24_stand, aes(x= D_oxygen_mgL, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1, alpha = 0.5)+
labs(x = "Dissolved oxygen (mgL)", y = "Crabs consumed") + Theme1 +
scale_x_continuous(breaks = c(0.1082296, 0.384854, 0.66147994, 0.935517247), labels = c("4", "6", "8", "10"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,8), breaks = seq(0, 6, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("OxC.jpeg", units="mm", width=90, height=60, dpi=300)
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
4) Bare rock Crabs 24 h
ggplot() +
geom_ribbon(data = p_C24BR, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_C24BR, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Crabs24_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
labs(x = "Bare rock cover (%)", y = "Consumed crabs") + Theme1 +
scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("BareC.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_point()`).
4) Fish length - Crabs 24 h
ggplot() +
geom_ribbon(data = p_C24FL, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_C24FL, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Crabs24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
geom_jitter(data = Crabs2_stand, aes(x= Bare_rock, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
labs(x = "Fish length (cm)", y = "Consumed crabs") + Theme1 +
scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("Fish_length.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 7 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
None
p_S24T <- ggpredict(my_glmS24_T1, terms ="Temperature [all]")
p_S24S <- ggpredict(my_glmS24_S1, terms ="Salinity_ppt [all]")
1) Temperature vs squidpops
ggplot() +
geom_ribbon(data = p_S24T , aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_S24T , aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Squidpops24_stand, aes(x= Temperature, y=Consumed_individuals), color = "orange1", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
geom_jitter(data = Squidpops2_stand, aes(x= Temperature, y=Consumed_individuals), color = "darkslateblue", width=0.02, height = 0.1, size = 1.5, alpha = 0.5)+
labs(x = "Temperature (°C)", y = "Consumed squidpops") + Theme1 +
scale_x_continuous(breaks = c(-0.04868, 0.21096, 0.4706, 0.73023, 0.98987), labels = c("10", "12", "14", "16", "18"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("TemperatureS.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
2) Salinity vs squidpops
ggplot() +
geom_ribbon(data = p_S24T, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_S24T, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Squidpops24_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "orange1", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
geom_jitter(data = Squidpops2_stand, aes(x= Salinity_ppt, y=Consumed_individuals), color = "darkslateblue", width=0.01, height = 0.1, size = 1.5, alpha = 0.5)+
labs(x = "Salinity (ppt)", y = "Consumed squidpops") + Theme1 +
scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
scale_y_continuous(limits = c(0,5.8), breaks = seq(0, 5, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 21 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("SalS.jpeg", units="mm", width=110, height=80, dpi=300)
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 30 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot() +
geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
labs(x = "Bare rock cover (%)", y = "Consumed squidpops") + Theme1 +
scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("BRS.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot() +
geom_jitter(data = Squidpops2_stand, aes(x= Fish_size, y=Consumed_individuals), color = "darkslateblue", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
geom_jitter(data = Squidpops24_stand, aes(x= Fish_size, y=Consumed_individuals), color = "orange1", width=0.015, height = 0.1, size = 1.5, alpha = 0.5)+
labs(x = "Fish Length (cm)", y = "Consumed squidpops") + Theme1 +
scale_x_continuous(breaks = c(0.011979959
,0.165667696
, 0.319355434
, 0.473043171
, 0.626730908
, 0.780418645
,0.934106383), labels = c("14","16", "18", "20", "22", "24", "26"), limits = c(0,1)) +
scale_y_continuous(limits = c(-0.5,6), breaks = seq(0, 5, 1)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 12 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("FLS2.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
p_Lphloro <- ggpredict(my_glmL_PL1, terms ="phloro_mean [all]")
p_L24S <- ggpredict(my_glmL_S1, terms ="Salinity_ppt [all]")
p_L24UA <- ggpredict(my_glmL_UA1, terms ="Understorey_algae [all]")
ggplot() +
geom_ribbon(data = p_Lphloro, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_Lphloro, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Herb_all_stand, aes(x= phloro_mean, y=Consumed_proportion), color = "orange1", width=0.02, height = 0.01, size = 1.5, alpha = 0.5)+
labs(x = "Phlorotannin content (mg"~g^-1*" DW)", y = "Blade consumed proportion") + Theme1 +
scale_x_continuous(breaks = c(0.04788192, 0.2840443, 0.5202066, 0.756369, 0.9925314), labels = c("8", "16", "24", "32", "40"), limits = c(0,1)) +
scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 230 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("PhloroL.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 229 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot() +
geom_ribbon(data = p_L24S , aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_L24S, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Herb_all_stand, aes(x= Salinity_ppt, y=Consumed_proportion), color = "orange1", width=0.02, height = 0.01, size = 1.5, alpha = 0.5)+
labs(x = "Salinity (ppt)", y = "Blade consumed proportion") + Theme1 +
scale_x_continuous(breaks = c(0.09803922, 0.49019608, 0.88235299, 1.2745098), labels = c("33", "33.5", "34", "33.5"), limits = c(0,1)) +
scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 268 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("SalL.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 240 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot() +
geom_ribbon(data = p_L24UA, aes(x = x, y = predicted, ymin = conf.low, ymax = conf.high), fill = "orange1", alpha = .25) +
geom_line(data = p_L24UA, aes(x = x, y = predicted), color = "orange1", size = 1) +
geom_jitter(data = Herb_all_stand, aes(x= Understorey_algae, y=Consumed_proportion), color = "orange1", width=0.03, height = 0.03, size = 1.5, alpha = 0.5)+
labs(x = "UA cover (%)", y = "Predicted consumed proportion ") + Theme1 +
scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1), labels = c("0", "20", "40", "60", "80", "100"), limits = c(0,1)) +
scale_y_continuous(limits = c(0,1.1), breaks = seq(0, 1, 0.2)) +
theme(legend.position = c(0.1, 0.85))
## Warning: Removed 269 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave("UAL.jpeg", units="mm", width=90, height=70, dpi=300)
## Warning: Removed 275 rows containing missing values or values outside the scale range
## (`geom_point()`).
Comparison between data obtained in-situ and long-term data. Sea surface temperature (SST) and dissolved oxygen daily average were extracted from E.U. Copernicus Marine Service Information (doi =10.48670/moi-00169, doi =10.48670/moi-00015, respectively). Salinity data were extracted using the NASA Earth Data Search (doi = SMP20-4U7CS). SST and salinity data were collected between 1 December 2018 and 31 March 2019, which matches the period of the study. Dissolved oxygen data were available for the time period 1 November 2021 and 28 February 2025.
Temp <- ggplot() +
geom_jitter(data = temp_data, aes(x= Latitude*-1, y=analysed_sst), width = 0.08, height = 0.03, size = 0.8, alpha = 0.3) +
geom_boxplot(data = temp_data, aes(x= Latitude*-1, y=analysed_sst, group = cut_width(Latitude,1)), alpha = 0.2) +
geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$Temperature), size=2, color = "red", pch =15)+
Theme1 +
labs( y="SST (°C)", x="Latitude (°S)") +
theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x = element_text(size = 11))
Temp
ggsave("TVar.jpeg", units="mm", width=130, height=100, dpi=300)
Sal <- ggplot() +
geom_jitter(data = Sal_data, aes(x= lat*-1, y=Salinity), width = 0.08, height = 0.03, size = 0.8, alpha = 0.6) +
geom_boxplot(data = Sal_data, aes(x= lat*-1, y=Salinity, group = cut_width(lat,1)), alpha = 0.2) +
geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$Salinity_ppt), size=2, color = "red", pch =15)+
Theme1 +
labs( y="Salinity (ppt)", x="Latitude (°S)") +
theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x = element_text(size = 11))
Sal
ggsave("SalVar.jpeg", units="mm", width=130, height=100, dpi=300)
library(respR)
## Warning: package 'respR' was built under R version 4.4.2
## Registered S3 method overwritten by 'respR':
## method from
## print.inspect pryr
##
## Attaching package: 'respR'
## The following object is masked from 'package:MASS':
##
## select
o2_data$conv <- convert_DO(o2_data$o2, # data to convert
from = "umol/L", # oxygen unit to convert from
to = "mg/L", # oxygen unit to convert to
t = 15, # in C
S = 34, # in ppt
P = 1.013) # in bar
ox <- ggplot() +
geom_jitter(data = o2_data, aes(x= Latitude*-1, y=conv), width = 0.08, height = 0.03, size = 0.8, alpha = 0.2) +
geom_boxplot(data = o2_data, aes(x= Latitude*-1, y=conv, group = cut_width(Latitude,1)), alpha = 0.2) +
geom_point(data =Environmental_sum, aes(x= Latitude*-1,y=Environmental_sum$D_oxygen_mgL), size=2, color = "red", pch =15)+
Theme1 +
labs( y="Dissolved oxygen (mg"~ L^-1*")", x="Latitude (°S)") +
theme(axis.text.y = element_text(size = 11), axis.title.y = element_text(size = 11), axis.text.x = element_text(size = 11), axis.title.x = element_text(size = 11))
ox
ggsave("o2Var.jpeg", units="mm", width=130, height=100, dpi=300)
In this section, a PCA with all the predictor variables was carried out and contracted with latitude.
Variables <- Phlorotannins_sum %>%
left_join(Environmental_sum, by = "Site") %>%
left_join(Group_abundances_mean, by = "Site") %>%
left_join(Richness_meanpred, by = "Site") %>%
left_join( pcoa_df1, by = "Site") %>%
left_join(Size_mean_wide, by = "Site") %>%
left_join(Benthic_sum_wide_var, by = "Site")
Variables1 <- dplyr::select(Variables, -Site, -Latitude, -Day_length)
# Perform PCA
pca_result <- prcomp(Variables1, center = TRUE, scale. = TRUE)
# View PCA summary
summary(pca_result)
## Importance of components:
## PC1 PC2 PC3 PC4 PC5 PC6 PC7
## Standard deviation 1.9056 1.8232 1.5595 1.2683 1.1188 0.80575 0.67715
## Proportion of Variance 0.2594 0.2374 0.1737 0.1149 0.0894 0.04637 0.03275
## Cumulative Proportion 0.2594 0.4968 0.6705 0.7854 0.8748 0.92121 0.95396
## PC8 PC9 PC10 PC11 PC12
## Standard deviation 0.59107 0.42934 0.32003 0.09201 4.97e-15
## Proportion of Variance 0.02495 0.01317 0.00732 0.00060 0.00e+00
## Cumulative Proportion 0.97891 0.99208 0.99940 1.00000 1.00e+00
# Access PCA components
pca_result$rotation # Loadings (eigenvectors)
## PC1 PC2 PC3 PC4
## phloro_mean 0.15312183 -0.38036071 0.16987612 -0.40416794
## D_oxygen_mgL 0.04416129 -0.06947408 0.22900240 -0.29578882
## Salinity_ppt 0.32770211 -0.28725268 -0.30041598 0.01204521
## Temperature 0.33959317 -0.19864160 0.25173653 0.07112323
## V_Visibility -0.26155336 -0.04655336 -0.46754924 -0.10336321
## Herbivores 0.17945120 0.26495508 -0.32071365 -0.30175664
## Invertebrate_predators 0.23631635 0.35158153 0.08199631 -0.22146569
## Predatory_fish 0.33348022 0.19557271 0.12844257 0.42742255
## Mean_SR 0.24588446 0.13891112 -0.46658761 -0.10235712
## PCo1 -0.45336621 -0.12321574 0.20140955 0.11357894
## Fish_size 0.40286210 -0.15837393 0.19130840 -0.19017485
## Kelp 0.20344908 -0.25734537 -0.07308674 0.49797138
## Bare_rock -0.02538404 0.41986093 0.34602767 -0.18266636
## Understorey_algae -0.11785619 -0.43843959 -0.02128394 -0.26379807
## PC5 PC6 PC7 PC8
## phloro_mean -0.22928655 -0.061934912 0.16030316 0.06977640
## D_oxygen_mgL 0.71917764 -0.141970506 0.24870627 0.01063205
## Salinity_ppt -0.04621847 0.088594787 -0.25384253 0.37791097
## Temperature 0.34233239 0.366339285 0.01239502 0.29709948
## V_Visibility 0.28277131 -0.265296643 -0.32346405 -0.03883120
## Herbivores -0.18909046 -0.075864110 0.65389541 0.11582022
## Invertebrate_predators 0.07850895 -0.516299347 -0.21984644 0.32372054
## Predatory_fish -0.04690808 -0.353793051 -0.09874146 0.21766508
## Mean_SR 0.31766340 0.122121426 -0.05131298 -0.29528597
## PCo1 0.18079042 -0.201314998 0.07870739 0.07955889
## Fish_size -0.07309108 -0.101363639 -0.28571931 -0.59602278
## Kelp 0.08064145 -0.390533765 0.39301185 -0.31521002
## Bare_rock -0.07793862 -0.008588723 -0.09287449 -0.20915738
## Understorey_algae -0.18785568 -0.384416932 -0.06340740 0.07047988
## PC9 PC10 PC11 PC12
## phloro_mean -0.19914246 -0.302226749 0.547906843 -0.325191416
## D_oxygen_mgL 0.28917754 -0.203944032 -0.080842398 -0.005764669
## Salinity_ppt 0.29382455 -0.354717840 0.055501156 0.535005519
## Temperature -0.10134383 0.454952128 -0.004250124 -0.076982913
## V_Visibility 0.07492272 -0.003118584 0.171242875 -0.276713426
## Herbivores 0.13103355 -0.015917972 -0.278008088 0.010513205
## Invertebrate_predators -0.52190028 0.112816781 0.002910929 0.218806222
## Predatory_fish 0.38672379 -0.153625509 -0.023726204 -0.490900382
## Mean_SR -0.06100320 0.112910544 0.207166319 -0.035055321
## PCo1 -0.11580691 -0.293793915 -0.143554874 0.226404458
## Fish_size -0.04289802 -0.178281745 -0.435374536 0.012026704
## Kelp -0.11866981 0.127147103 0.290625846 0.278306208
## Bare_rock 0.44193909 0.139799385 0.473690956 0.325001830
## Understorey_algae 0.32715632 0.575421427 -0.145109455 0.028404422
pca_result$x # Principal components (scores)
## PC1 PC2 PC3 PC4 PC5 PC6
## [1,] 1.01934232 -2.465942186 -0.05856192 -1.1805419 -0.74403145 -0.6008792
## [2,] 1.64364575 -0.247652503 0.33059729 2.8574056 0.41942831 -0.9032240
## [3,] 0.65673090 -1.475702438 -0.83220277 0.1380973 -0.43291297 1.0595816
## [4,] 0.10429367 3.053015250 -2.55270043 -1.8166224 0.04174398 -0.5924529
## [5,] 0.99165334 0.069541318 0.37106187 -0.4478629 -1.35364440 0.5912784
## [6,] -4.14190953 -0.901118939 -2.06185691 1.0229257 -0.76843221 -0.5500265
## [7,] 1.57351825 1.793784930 1.24878295 -0.1266223 -0.94011881 0.5183558
## [8,] 0.52079420 2.322458024 1.45717978 0.2800948 0.28004323 -1.0837141
## [9,] -0.02668262 -0.009488045 -0.64428859 -0.2340490 2.98842103 0.7665162
## [10,] -0.39575077 -2.747177050 1.78174925 -1.6003581 0.71592698 -0.7627492
## [11,] -3.50506614 1.243557346 2.39160564 0.2027377 -0.10405266 0.8614225
## [12,] 1.55943063 -0.635275706 -1.43136615 0.9047954 -0.10237103 0.6958914
## PC7 PC8 PC9 PC10 PC11 PC12
## [1,] 0.2498895 -0.72786447 0.49737574 0.47748733 0.08635639 -4.642120e-15
## [2,] 0.5267490 -0.03407682 -0.56423696 0.21510138 0.01157954 -5.436623e-15
## [3,] -0.8074396 0.10588014 -0.19258732 0.36593723 -0.18791464 -5.256212e-15
## [4,] 0.6350624 -0.05283392 -0.36977732 0.09519209 -0.04917014 -4.992534e-15
## [5,] -0.6701076 -0.56860844 -0.67721358 -0.40815778 0.11667248 -4.359360e-15
## [6,] -0.4861734 0.48875365 0.14721247 -0.08360601 0.05770054 -5.280498e-15
## [7,] 0.1000942 1.32580148 0.29863352 0.19812510 0.07925996 -4.798245e-15
## [8,] -1.1046837 -0.41762546 0.51534582 -0.17420288 -0.07669713 -4.576201e-15
## [9,] -0.4462995 -0.02123266 0.04596192 0.08061882 0.10601955 -4.909267e-15
## [10,] 0.3077744 0.66291487 -0.28732577 -0.34383257 -0.06506740 -4.635181e-15
## [11,] 0.8137796 -0.57304152 -0.01998176 0.14372574 -0.03156541 -4.538037e-15
## [12,] 0.8813548 -0.18806684 0.60659323 -0.56638845 -0.04717373 -4.943962e-15
# Scree plot
fviz_eig(pca_result, addlabels = TRUE)
pca_scores <- as.data.frame(pca_result$x)
pca_scores$Site <- Variables$Site
pca_scores$Latitude <- Variables$Latitude
pca_loadings <- as.data.frame(pca_result$rotation)
pca_loadings$Variable <- rownames(pca_loadings)
# Rename variable labels in the pca_loadings data frame
pca_loadings$Variable <- dplyr::recode(pca_loadings$Variable,
V_Visibility = "Visibility",
Understorey_algae = "UA cover",
Bare_rock = "Bare rock cover",
D_oxygen_mgL = "Dissolved oxygen",
phloro_mean = "Phlorotannins",
Kelp = "Kelp cover",
PCo1 = "PCo1",
Temperature = "Temperature",
Fish_size = "Fish length",
Salinity_ppt = "Salinity",
Predatory_fish = "Predatory fish",
Invertebrate_predators = "Invertebrate predators",
Mean_SR = "Species richness",
Herbivores = "Herbivores"
)
ggplot() +
# Add PCA scores (samples) with latitude gradient
geom_point(data = pca_scores, aes(x = PC1, y = PC2, color = Latitude*-1), size = 3) +
scale_color_viridis_c(option = "plasma", trans = "reverse") +
# Add PCA loadings (variables as arrows)
geom_segment(data = pca_loadings,
aes(x = 0, y = 0, xend = PC1 * 5, yend = PC2 * 5),
arrow = arrow(length = unit(0.2, "cm")), color = "black") +
# Add variable labels
geom_text(data = pca_loadings,
aes(x = PC1 * 5, y = PC2 * 5, label = Variable),
color = "black", vjust = -0.5) +
# Customize plot appearance
labs( x = "Principal Component 1",
y = "Principal Component 2",
color = "Latitude (°S)") +
Theme1
ggsave("PCA.jpeg", units="mm", width=130, height=120, dpi=300)
# Load vegan package
library(vegan)
# Create a distance matrix from PCA scores
pca_dist <- dist(pca_scores[, c("PC1", "PC2")])
# Perform PERMANOVA
permanova_result <- adonis2(pca_dist ~ Latitude, data = pca_scores, permutations = 999)
print(permanova_result)
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = pca_dist ~ Latitude, data = pca_scores, permutations = 999)
## Df SumOfSqs R2 F Pr(>F)
## Latitude 1 33.467 0.43742 7.7754 0.001 ***
## Residual 10 43.042 0.56258
## Total 11 76.510 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A summary table of the size estimations by species and site was created.
Size_summary <- Size_metadata %>%
group_by(Species) %>%
summarise (N_individuals = n(), mean_Size_sp = mean(Size_cm), sd = sd(Size_cm), Max_size = max(Size_cm), Min_size = min(Size_cm))
kable(Size_summary)
| Species | N_individuals | mean_Size_sp | sd | Max_size | Min_size |
|---|---|---|---|---|---|
| Acanthistius_pictus | 3 | 16.000000 | 1.0000000 | 17 | 15 |
| Anisotremus_scapularis | 4 | 21.250000 | 2.9860788 | 25 | 18 |
| Aplodactylus_punctatus | 54 | 27.944444 | 6.5196848 | 40 | 14 |
| Auchenionchus_critinus | 1 | 12.000000 | NA | 12 | 12 |
| Auchenionchus_sp. | 14 | 19.071429 | 6.0949448 | 32 | 11 |
| Auchenionchus_variolosus | 3 | 28.666667 | 5.6862407 | 35 | 24 |
| Bovichtus_chilensis | 3 | 17.666667 | 9.0184995 | 27 | 9 |
| Calliclinus_geniguttatus | 19 | 10.157895 | 2.3632172 | 14 | 4 |
| Cheilodactylus_variegatus | 125 | 22.560000 | 7.0715240 | 38 | 6 |
| Cheilotrema_fasciatum | 1 | 28.000000 | NA | 28 | 28 |
| Chromis_crusma | 85 | 13.611765 | 4.2457096 | 26 | 5 |
| Congiopodus_peruvianus | 2 | 13.500000 | 0.7071068 | 14 | 13 |
| Girella_laevifrons | 11 | 26.272727 | 7.3087743 | 38 | 15 |
| Graus_nigra | 15 | 25.066667 | 4.9923751 | 37 | 19 |
| Hemilutjanus_macrophthalmos | 2 | 20.000000 | 2.8284271 | 22 | 18 |
| Hepatus_chilensis | 1 | 6.000000 | NA | 6 | 6 |
| Homalaspis_plana | 20 | 9.800000 | 1.9084301 | 12 | 4 |
| Hypsoblennius_sordidus | 5 | 7.800000 | 1.0954451 | 9 | 6 |
| Isacia_conceptionis | 7 | 25.571429 | 2.5071327 | 30 | 22 |
| Labrisomus_philippii | 18 | 28.055556 | 6.9745101 | 38 | 15 |
| Nexilosus_latifrons | 4 | 23.000000 | 5.7154761 | 30 | 17 |
| Paralabrax_humeralis | 14 | 26.571429 | 8.4826908 | 44 | 15 |
| Paralichthys_microps | 4 | 32.500000 | 23.1012265 | 65 | 16 |
| Paraxanthus_barbiger | 12 | 6.250000 | 1.6583124 | 9 | 4 |
| Patagonotothen_brevicauda | 10 | 8.100000 | 1.3703203 | 10 | 6 |
| Pinguipes_chilensis | 81 | 28.370370 | 8.5197483 | 51 | 13 |
| Prolatilus_jugularis | 20 | 24.000000 | 6.9357958 | 39 | 10 |
| Romaleon_setosum | 17 | 8.058824 | 2.2768012 | 12 | 5 |
| Scartichthys_spp. | 176 | 13.670454 | 5.0927884 | 34 | 4 |
| Schroederichthys_chilensis | 3 | 51.666667 | 2.8867513 | 55 | 50 |
| Sebastes_oculatus | 12 | 12.333333 | 7.1137937 | 24 | 4 |
| Semicossyphus_darwini | 6 | 26.500000 | 6.4420494 | 36 | 20 |
| Taliepus_sp. | 27 | 9.444444 | 1.8045526 | 12 | 5 |
write.csv(Size_summary, "Size_summary2.csv")
In this part, information about the operating system and R packages attached during the production of this document can be found
sessionInfo() %>% pander
R version 4.4.1 (2024-06-14 ucrt)
Platform: x86_64-w64-mingw32/x64
locale: LC_COLLATE=English_Australia.utf8, LC_CTYPE=English_Australia.utf8, LC_MONETARY=English_Australia.utf8, LC_NUMERIC=C and LC_TIME=English_Australia.utf8
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: respR(v.2.3.3), fitdistrplus(v.1.2-2), survival(v.3.6-4), factoextra(v.1.0.7), lme4(v.1.1-35.4), Matrix(v.1.7-0), glmmTMB(v.1.1.10), mgcv(v.1.9-1), nlme(v.3.1-167), geosphere(v.1.5-18), MuMIn(v.1.48.4), MASS(v.7.3-60.2), corrplot(v.0.95), car(v.3.1-2), carData(v.3.0-5), knitr(v.1.49), performance(v.0.13.0), ggeffects(v.1.7.0), pander(v.0.6.5), patchwork(v.1.3.0), vegan(v.2.6-6.1), lattice(v.0.22-6), permute(v.0.9-7), summarytools(v.1.0.1), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.4), purrr(v.1.0.2), readr(v.2.1.5), tidyr(v.1.3.1), tibble(v.3.2.1), ggplot2(v.3.5.1) and tidyverse(v.2.0.0)
loaded via a namespace (and not attached): rstudioapi(v.0.16.0), DHARMa(v.0.4.6), jsonlite(v.1.9.0), shape(v.1.4.6.1), datawizard(v.1.0.0), magrittr(v.2.0.3), TH.data(v.1.1-2), estimability(v.1.5.1), magick(v.2.8.3), farver(v.2.1.2), nloptr(v.2.1.0), rmarkdown(v.2.29), ragg(v.1.3.2), vctrs(v.0.6.5), minqa(v.1.2.7), base64enc(v.0.1-3), rstatix(v.0.7.2), htmltools(v.0.5.8.1), haven(v.2.5.4), broom(v.1.0.6), marelac(v.2.1.11), sass(v.0.4.9), bslib(v.0.9.0), plyr(v.1.8.9), sandwich(v.3.1-1), emmeans(v.1.10.2), zoo(v.1.8-12), cachem(v.1.1.0), TMB(v.1.9.16), lifecycle(v.1.0.4), pkgconfig(v.2.0.3), sjlabelled(v.1.2.0), R6(v.2.6.1), fastmap(v.1.2.0), rbibutils(v.2.2.16), digest(v.0.6.35), numDeriv(v.2016.8-1.1), colorspace(v.2.1-0), textshaping(v.0.4.0), ggpubr(v.0.6.0), labeling(v.0.4.3), timechange(v.0.3.0), abind(v.1.4-5), compiler(v.4.4.1), withr(v.3.0.2), backports(v.1.5.0), ggsignif(v.0.6.4), gsw(v.1.2-0), roll(v.1.1.7), tools(v.4.4.1), glue(v.1.7.0), grid(v.4.4.1), checkmate(v.2.3.1), cluster(v.2.1.6), reshape2(v.1.4.4), see(v.0.10.0), generics(v.0.1.3), gtable(v.0.3.6), tzdb(v.0.4.0), data.table(v.1.15.4), hms(v.1.1.3), sp(v.2.1-4), xml2(v.1.3.6), ggrepel(v.0.9.6), pillar(v.1.10.1), splines(v.4.4.1), pryr(v.0.1.6), oce(v.1.8-3), tidyselect(v.1.2.1), reformulas(v.0.4.0), stats4(v.4.4.1), xfun(v.0.51), rapportools(v.1.1), matrixStats(v.1.5.0), stringi(v.1.8.4), yaml(v.2.3.8), boot(v.1.3-30), evaluate(v.1.0.3), codetools(v.0.2-20), tcltk(v.4.4.1), seacarb(v.3.3.3), SolveSAPHE(v.2.1.0), cli(v.3.6.2), RcppParallel(v.5.1.10), xtable(v.1.8-4), systemfonts(v.1.1.0), Rdpack(v.2.6), munsell(v.0.5.1), jquerylib(v.0.1.4), Rcpp(v.1.0.12), coda(v.0.19-4.1), parallel(v.4.4.1), bayestestR(v.0.15.1), viridisLite(v.0.4.2), mvtnorm(v.1.2-5), scales(v.1.3.0), insight(v.1.0.2), crayon(v.1.5.3), rlang(v.1.1.4) and multcomp(v.1.4-26)